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ABSTRACT 


A reduced shallow water model under constant, non-zero advection in infinite do¬ 
mains is considered. High-Order Givoli-Neta (G-N) and Hagstrom-Hariharan (H-H) non¬ 
reflecting boundary conditions (NRBCs) are introduced to create a finite computational 
space and solved using a spectral element formulation with high-order time integration. 
Numerical examples are used to demonstrate the synergy of using high-order spatial, time 
and boundary discretizations. Several alternatives are also presented for solving open do¬ 
main problems. These alternatives include adjustments to the G-N NRBC based on phys¬ 
ical arguments as well as formulating the boundary condition for arbitrary domains using 
unstructured grids. The H-H polar NRBC is also formulated in an unstructured grid setting 
and extended to include dispersive effects. Results show that by balancing all numerical 
errors involved, high-order accuracy can be achieved for unbounded channel problems. 
Further, the adjustments to the G-N and H-H NRBCs to operate in an unstructured grid 
setting are shown to significantly reduce errors over first order non-reflecting boundary 
schemes when operating in an open domain configuration. 
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I. INTRODUCTION 


Simulation of wave propagation through large - perhaps unbounded - domains has 
been an active area of research for several decades. Such studies are important to many 
applications such as acoustics, electromagnetics, meteorology, solid geophysics and aero¬ 
dynamics to name just a few. The combination of the complexity of the partial differential 
equation sets involved and the infinite possibilities of initial data mandates the numerical 
solution to such problems. Of course, to undertake the numerical solution on an infinite 
domain would be both foolhardy and impossible. Generally speaking, to overcome this 
computational challenge, it is quite common to truncate the infinite domain by imposing 
some type of boundary condition on a “sufficiently large” sub-domain that captures the area 
of interest. 

When truncating the domain, the modeler must devise boundary conditions for the 
truncated domain. Of course, by imposing a boundary where one does not physically 
exist, the problem is changed - and unless chosen carefully, would certainly be expected to 
pollute the solution as the problem evolves and impinges on the boundary. Therefore, two 
main possibilities exist for the modeler: 

• Choose a convenient, easily implementable boundary condition that does not nec¬ 
essarily reflect the physical problem and solve it on a large sub-domain. The idea 
behind this technique is that the boundary effects are negligible for a short time evo¬ 
lution of the problem in a small area of interest away from the boundaries. 

• Choose a boundary condition that preserves the true behavior of the infinite solution 
at the boundary and solve the problem on a smaller sub-domain. The idea behind this 
technique is that the additional effort extended to impose a better boundary condition 
will be worth the effort and allow for solving the problem on a smaller domain. 

For obvious reasons, the first possibility has only limited usefulness. To see why, suppose 
that we wanted to model the wave motion following a pebble dropped in the center of a 
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large, still pond. Now, suppose that we have a truncated domain to model this phenomena 
- say a bathtub. If the pebble is dropped in the bathtub, the waves generated by the pebble 
would propagate much like that in the pond - until, that is - the wave front reaches the hard 
walls of the bathtub. At this point, the bathtub model ceases to be a useful representation 
of the pond due to behavior caused by the non-physical boundary. If the modeler wishes to 
see what happens a short time later - a larger bathtub would be required. This same prin¬ 
ciple would apply for the numerical solution of this propagation problem - a poor choice 
of boundary condition mandates the use of a larger computational domain. This in turn, 
requires additional computational resources. For this reason, much effort has and continues 
to be exerted on finding suitable boundary conditions that apply on smaller domains. 

This dissertation examines the use of high-order non-reflecting boundary conditions 
(NRBCs) to solve a class of infinite domain, wave propagation problems. In the last 35 
years or so, much research has been done to develop NRBCs that, after discretization, lead 
to a scheme that is stable, accurate, efficient and easy to implement. Of course, it is difficult 
to find a single NRBC that is ideal in all respects and all cases; this is why the quest for 
better NRBCs and their associated discretization schemes continues. 

Sequences of increasing-order NRBCs have been available for a long time (e.g., the 
Bayliss-Turkel conditions [1] constitute such a sequence), but they had been regarded as 
impractical beyond 2nd or 3rd order from the implementation point of view. Only since the 
mid 90s have practical high-order NRBCs been devised. The first high-order local NRBC 
was proposed by Collino [2], for two-dimensional time-dependent waves in rectangular 
domains. Its construction requires the solution of the one-dimensional wave equation on 
the boundary. Grote and Keller [3] developed a high-order converging NRBC for the three- 
dimensional time-dependent wave equation, based on spherical harmonic transformations. 
Sofronov [4, 5] proposed exact boundary conditions for the three- and two- dimensional 
wave equations in spherical and polar coordinates, respectively (it is proved that NRBCs 
demonstrated in [3] and [4] are reduced to each other). 
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Hagstrom and Hariharan [6] constructed high-order NRBCs for the two- and three- 
dimensional time-dependent wave equations based on the analytic series representation for 
the outgoing solutions of these equations. For time-dependent waves in a two-dimensional 
wave guide, Guddati and Tassoulas [7] devised a high-order NRBC by using rational ap¬ 
proximations and recursive continued fractions. Givoli [8] has shown how to derive high- 
order NRBCs for a general class of wave problems, leading to a symmetric FE formulation. 
In [9], this methodology was applied to the particular case of time-harmonic waves, using 
optimally localized Dirichlet-to-Neumann (DtN) NRBCs (see also [10]). 

Previous studies of this nature have encountered limits in the accuracy of such solu¬ 
tions. These accuracy limits can be caused by time and space discretization as well as from 
the boundary scheme used in the solution of the problem. We seek to employ high-order 
numerical methods in time and space to diminish the effects of discretization error in order 
to determine the true efficacy of a given non-reflecting boundary condition. 

The rest of the dissertation is outlined as follows. Chapter II motivates and derives 
the equations for the problem under consideration. In Chapter III, we summarize the main 
boundary schemes currently in use and specifically show the Givoli-Neta (G-N) NRBC 
in detail. We then describe the high-order spectral element method used to discretize the 
problem in space (up to 16 th order polynomials) in Chapter IV. Chapter V discusses the 
Runge-Kutta time discretization demonstrated up to 10 t/l order. Chapter VI provides nu¬ 
merical examples in various configurations and conditions to demonstrate concepts. Chap¬ 
ter VII considers challenges associated with arbitrary boundary configurations and provides 
results for a low-order boundary treatment using unstructured grids. Chapter VIII exam¬ 
ines the potential for exploiting other non-reflecting boundary conditions employed in an 
unstructured grid formulation and concludes with areas of future research. 
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II. EQUATIONS OF FLUID MOTION 


In order to scope the enormous problem of wave propagation and provide a good 
test bed of examples for simulation of the non-reflecting boundary conditions under exami¬ 
nation, we will derive the shallow water equations. The set of equations under consideration 
have been used to predict Tsunamis and storm surges [11], as well as modeling atmospheric 
flows. The term “shallow water” is a bit deceiving, as the medium does not necessarily have 
to be water, nor does it have to be shallow. To elaborate - the equations under consideration 
are generated by very general physical principles, namely conservation of mass and mo¬ 
mentum, that are then simplified using reasonable assumptions. In this derivation, we will 
abbreviate the approach taken by Cushman-Roisin [12], Pedlosky [13] and Batchelor [14] 
and further simplify the equations by reducing them to a scalar Klein-Gordon equation 
equivalent. 

A. CONSERVATION OF MASS 

One of the fundamental physical principles is that mass can neither be created nor 
destroyed. Consider a control volume - a fixed region in space where fluid is allowed 
to occupy and pass through. Within this control volume, mass is conserved. In other 
words, for the mass to change in a control volume, there must be mass passing through the 
boundary of the control volume. Suppose dm is an infinitesimal portion of the mass, dV 
and dA are infinitesimal portions of the control volume and its boundary, respectively, and 
p is the density of the fluid occupying the control volume. Then by this argument, the mass 
enclosed by the surface at any instant is f dm — f p dV. Further, the net rate that mass 
is flowing outwards across the boundary is f (pu) ■ n dA. Putting these two ideas together 
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we find that for mass to be conserved: 


^ J pdV = - J (pu) • n dA. 

' s V ^ ^ V* ^ 

Time rate of mass Net rate of mass flux 

change in CV across boundary 

Here: u = (u. v. w) T is the fluid velocity and n is the outward pointing unit normal on the 
boundary. 

Further, upon differentiation under the integral sign (remembering that the control 
volume is fixed in space) and transforming the surface integral using the divergence theo¬ 
rem: 


/ ^ dV + j V ■ (pu) dV = 0 
/(| +v . pu ) dV=0 

This relation is valid for all choices of the control volume, therefore, the integrand 
must also be identically zero, i.e., 

^ + V • pu = 0 (III) 

The differential equation (II. 1) is commonly referred to as the continuity equation in fluid 
mechanics. 

B. IMPORTANCE OF THE EARTH’S ROTATION 

Modeling phenomena on a large scale, such as the currents of the ocean or winds in 
the atmosphere, may require special handling since the earth is not static. In fact, the earth 
rotates on its axis approximately once every 24 hours. This results in a mean rotation rate 
fl of: 

n = _^_. 

rotation period 
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Figure 1: Fixed (. X , Y) and rotating (x, y ) frameworks of reference. 

When the “exact” 1 rotation period of the earth is considered, this results in a mean angular 
velocity of 7.292 x 10~ 5 ra ^ ns [15]. The trajectory of a fluid in motion is expected to be 
influenced by this rotation if the fluid traveling at the speed U covers a distance L in a time 
interval greater than the rotation period. This concept is captured by a non-dimensional 
parameter £ and is defined as: 

27t/0 27 tU 

£ = L/U = ~VtL 

If e is on the order of or less than unity (e < 1), then we would expect that rotation 
is important. This number (neglecting the constant multiple 2n) is known as the Rossby 
number [12]. 

1. Equations in a Rotating Frame 

Now, we examine this rotating frame of reference on the earth. From the human 
perspective on the surface of the earth, we appear to be on a 2-dimensional surface. Suppose 
that we have X— and Y — axes that are the fixed or inertial reference frame and x — and 
y— axes that form the same reference frame, but rotating at the angular rate of Q. The unit 
vectors that follow this convention are defined by (I,J) and (i,j) as shown in Figure l 2 . It 

'The mean day is 23 hours, 56 minutes and 4.091 seconds, but variations caused by friction from the 
earth’s tides, as well as significant geophysical events on earth have been observed to cause fluctuations in 
this measure. 

2 Adapted from [12], Figure 2-1, p. 17. 
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follows that 


i = I cos fit + J sin f It j = —I sin fit + J cos fit 

and the coordinates of the position vector r = XI + Y J = xi + yj are correspondingly 

x = X cos fit + Y sin fit 
y = —X sin fit + Y cos fit 

Differentiating once with respect to time gives the rate of change of the coordinates relative 
to the moving frame, u = ^ = ui + nj (relative velocity). Differentiating again with 
respect to time gives the rate of change of the relative velocity in the moving frame, a = 
= ai + bj (relative acceleration). When completed and simplified, the absolute 
acceleration in the inertial frame with respect to the relative acceleration is: 

A = (a — 2fln — fl 2 a;) i + (b + 2flw — fl 2 r/) j (II.2) 

= Ai + Bj 

These results could also be derived in a vector form as outlined in [13] by defining 
the vector rotation in a direction common to both the rotating and inertial frames of refer¬ 
ence - f2 = flk, where k is a unit vector orthogonal to the plane. In this case, we can write 
(II.2) as: 

A = a + 2flxu + flx(flxx) (II.3) 

where u is extended to u = wi-Hy +tck. It should be noted that there are three contributions 
to the acceleration in the rotating frame: relative acceleration (a), one proportional to fl and 
the velocity, and one proportional to fl 2 and the position. The contribution proportional to 
fl and the velocity is known as the Coriolis acceleration and the other, proportional to Q 2 
and the position is the centripetal acceleration. 
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For practical purposes, the centripetal acceleration terms are often neglected since 
fl 2 ~ O (1CT 9 ). Additionally, even though centripetal acceleration causes objects on the 
surface of the planet to feel an outward pull, these objects do not fly off into space. In 
fact, even objects at rest, (u, v ) = 0 thus removing the Coriolis effect, for all intents and 
purposes, remain at rest as the gravitational pull of the earth keeps centripetal acceleration 
in check. 

When we neglect the centripetal acceleration terms in (II.3), the absolute accelera¬ 
tion terms in the inertial frame simplify to 

A = a + 2 fl x u 

= (a — 2Qv) i + (6 + 2fht) j + ck (II.4) 

2. 3-Dimensional Rotating Earth Model 

Now, we consolidate our results to apply to a 3-dimensional rotating earth model. 
Consider Figure 2 3 where U is oriented along the axis of rotation and an object is located 
on the surface at a latitude qb. A local coordinate system is set up with the axis orientation 
(x, y, z ) —» (east,north,radial) with standard convention of unit normal vectors. In this 
frame of reference, the earth’s rotation vector is expressed as 

f2 = cos 0j + fl sin 0k. (II.5) 

Using this to expand (II.4), we find that the acceleration in the inertial reference in 
terms of the rotating components has the following components: 


i : 

a + 2Qvj cos 0 — 2Qv sin 0 


j : 

b + 2f 2u sin 0 

(H.6) 

k : 

c — 2 Qu cos 0. 



3 Adapted from [12], Figure 2-8, p. 27. 
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Figure 2: Local Cartesian framework on spherical earth. 

Here, we notice the terms dependent on the latitude have common components, namely: 

/ = 2f2 sin (j) (II.7) 

/* = 2fi cos (f) 

The coefficient / is called the Coriolis parameter and the latter is the reciprocal Coriolis 
parameter. If we examine the scales of the components described in (II.6), we note that 
if we are describing geophysical fluid motion, the depth is much smaller than that of the 
surface it covers, and as such, the flows in this context tend to be parallel to the surface 
minimizing the effects of vertical flows. For our model, this implies that the effects of w 
are negligible compared to the effects of u and v. Additionally, any acceleration induced by 
the rotation in the vertical direction will be negligible compared to those along the surface. 
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This simplifies our acceleration in the inertial reference to: 


i : a — 2ttv sin 0 

j : b + 2 £lu sin 0 

k : c. 

In vector form, this can be written as: 

A = a + / (k x u) (II.8) 


Cushman-Rosin [12] provides two anecdotal justifications for this simplification. While 
the atmospheric layer that determines our weather is only about 10 km thick, cyclones and 
anticyclones spread over thousands of ki lometers. The second in the context of oceanic 
currents notes that flows are generally confined to the upper hundred meters but spread 
over tens of kilometers. 

C. CONSERVATION OF MOMENTUM 


The linear momentum of an object of mass m moving with velocity u is defined 
to be the product of the mass and velocity: P = mu, and in a closed system, must be 
conserved. Linear momentum is related to the forces acting using Newton’s second law of 
motion - namely, that the time rate of change of the momentum of an object is equal to the 
resultant forces on the object. That is, 


F = 



(H-9) 


This implies that if resultant forces are zero, the momentum of the particle must be constant. 
For this to happen, all of the forces (body and surface) acting on an object must sum to zero. 
In terms of the control volume framework discussed previously, this further implies that if 
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momentum is entering or leaving the control volume: 


d_ 

Ft 


p 


d 

= Ft mu * 


j pu*dV = I pbdV — I (pu)u-noL4 + f TncL4. 


Body forces 
acting on Q 


Net momentum flux 
across bdry of fl 


Forces acting 
on bdry of 


(II. 10) 


Here, u* is the velocity in the inertial frame. Body forces are those acting on the fluid 
volume that are proportional to the mass. The body forces considered here are gravity and 
(indirectly) the Coriolis force described in II. B. Others could include electromagnetic and 
centrifugal forces pertinent in alternate applications. 


1. Gravity Effects 

Gravity acts on a control volume strictly towards the center of the earth, and in the 
local coordinate system is along —k. This implies: 


pb g = -peg k 


(n.ii) 


Here, g is the gravitational constant which varies based on the distance from the center of 
the earth. At sea level, this value is approximately 9.798(5, an d in this context can be taken 
to be constant. 


2. Coriolis Effects 

Coriolis acts in a way to “adjust” the control volume’s acceleration when going 
from a rotating to an inertial frame. This term enters via the expression on the left hand 
side in the time change of momentum. Explicitly, this is 


Wt f puJV = 


dp <9u* 

dt^ + p -dt 


dV = 


^u^ + pA^j dV 

+ pa + pf (k x u) 


dV 


d 


= — / pudV + / p/(kxu) dV 
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provided that p is constant in time. 


3. Surface Force Effects 

Surface forces are those exerted across the boundary by the surrounding matter. 
Typical surface forces include pressure and viscosity. Since we are examining the equations 
in the context of atmospheric and oceanic applications, the effects of viscosity are small 
in comparison to other forces, and as such, will not be considered here. VanJoolen [16] 
derives these terms in detail for the interested reader, but further concludes the contribution 
of viscosity can be neglected in this application. The total boundary force exerted by the 
surrounding matter at the surface of a control volume is 

/ TneL4 = — [ pndA 


where p is the pressure exerted on the control volume by the surroundings. This surface in¬ 
tegral may be transformed to an integral over the volume by the analogue of the divergence 
theorem for a scalar quantity [14, p. 15] giving 

— J pndA = — J VpdV. 


Alternate derivations of this quantity can be found in [16] and [17], and provide further 
insights as from where this quantity is derived. 


4. Momentum Equations in a Rotating Frame 

We now apply the divergence theorem to transform the surface integral in (II. 10) 
and collect results. 


d_ 

dt 


pudV + 


V ■ (puu) dV 


pg\tdV + 


pf (k x u) dV + 


VpdV = 0 
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Here, uu is a tensor product of the velocity vector. Again, differentiating under the integral 
sign and consolidating terms yields: 


d . 

m {pu) 


V • (puu) — pgk + pf (k x u) + Vp 


dV = 0. 


Since this relation is valid for all choices of the control volume, the integrand must identi¬ 
cally be zero, i.e., 


— (pu) + V ■ (puu) + Vp = pgk + pf (u x k). 


( 11 . 12 ) 


Together with the continuity equation (II. 1) we have our equations of fluid motion in a 
rotating frame. 

D. FURTHER SIMPLIFICATIONS OF THE SYSTEM 


Our equations of fluid motion have taken several physical principles into account, 
but we wish to simplify them more in order to get a tractable model that can serve as a 
launching point for more testing. In the case of this analysis, we wish to make the following 
assumptions about the physical problem in order to arrive at the shallow water equations: 

The fluid is homogeneous: We assume that the fluid density is constant and uniform through¬ 
out the domain. 

The fluid is inviscid: This implies that the only surface force acting on the fluid is pressure 
(neglects shear forces which would act to retard the motion of the fluid.) 

The fluid is incompressible: Together with the assumption of homogeneity this decou¬ 
ples the dynamics from any thermodynamic considerations that might be used in 
another setting. 

Centrifugal forces are balanced by gravity: This allows simplification of acceleration terms 
in the rotating frame of reference. 
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Figure 3: The shallow water model with irregular bottom topography. 

The fluid is shallow: The depth terms in the applications considered are much smaller 
than the surface it covers, and therefore implies that flows are primarily along the 
surface. 

We consider a sheet of fluid as shown in Figure 3 with properties as outlined above. 
Here, we define the irregular bottom height below a sensible a reference value 2 = 0 as 
This reference level could be considered the fluid height at rest. The height of 
the surface of the fluid above the same reference level is defined as h(x, y, t). The depth 
of the fluid is therefore H(x, y, t) = h(x, y, t) + h B (x, y). To reiterate the shallow water 
assumption, we have a value for the height of the fluid II (x. y, t) that is much smaller than 
the length and width of the fluid. Again, we consider a fluid that is in the presence of 
rotation O about the z—axis. 

1. Mathematical Simplifications 

As outlined in [17, Appendix A], the momentum equations (11.12) can be mathe¬ 
matically simplified using nothing but product rule expansions and the continuity equation 
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to yield: 


du du du du 

ai + u dx + v dy 
dv dv dv 
di +u dx +v dy 
dw dw dw dw 

m +u ai + v d v +w al 


W 


w 


dz 

dv 

dz 


1 dp 

pdx~ V 
1 dp 

pdy = - fU 
1 dp__ 
p dz ^ 


(11.13) 


2. Implication of Homogeneity 

Since the fluid is assumed to be homogeneous in nature (p is constant), the conti¬ 
nuity equation (II. 1) reduces to 


du dv dw 
dx + dy + dz 


V ■ u = 0. 


(H.14) 


3. Implication of Shallow Fluid 

This assumption allows significant simplification of our fluid motion model. Here, 
we assume that the surface scale of the problem at hand is much larger than that of any 
depth considerations. Pedlosky [13] outlines a scaling argument that shows how the rel¬ 
ative importance of terms in the z—momentum equation allow all of the terms except for 
the pressure derivative and the gravity to be neglected. This collapses the z—momentum 
equation to 


dp 

dz 


-pg- 


We can then immediately depth integrate to yield 


P = ~pgz + A(x,y,t). 
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If the surface is under some constant ambient pressure p(x, y, h ) = p 0 , this implies that 
A(x, y, t) — p 0 + pgh(x, y, t) thereby, giving us an expression for p 


P(x, y, h) = pg (h{x, y, t) - z) + p 0 ■ 


We compute the pressure gradients in the x— and y— directions 


dp dh(x,y,t) dp dh(x,y,t ) 

dx ^ dx dy ^ dy 


(11.15) 


Here, we note [13] that the pressure gradients are independent of z so that the horizontal 
accelerations must also be independent of z. For consistency, we therefore assume that the 
horizontal velocities will also be independent of z. 


a. Primarily Horizontal Flow 

We have already observed that since the flow in shallow waters is primarily 
along the surface, the z— momentum collapsed down significantly. We can use this argu¬ 
ment to similarly simplify the x— and y— momentum equations. In this case, since w is 
significantly smaller than u or v. the terms w^ and w can also safely be neglected. Sub¬ 
stituting these simplifications along with our pressure gradients (11.15) into (11.13) results 
in the x— and y— momentum equations: 


du du du dh(x, y, t ) 

m + u dx + v dy - fv = 

dv dv dv dh(x, y, t ) 

di + U 9x + % ~ fU = — 


(11.16) 
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b. Continuity Equation Simplifications 

We now depth integrate our continuity equation (11.14) from z = — hb(x , y ) 
to z — h(x, y, t). After dropping variable dependencies for clarity, this yields: 


0 


V ■ u dz 


J — hs 

f (t + t) 

J- hB \dx dyj 


dz + w\ z=h - w\ z =- hB . 


(11.17) 


Considering reasonable boundary conditions for the last term as detailed in [12], we specify 
no normal flow on the rigid bottom (—hs) and a corresponding kinematic condition at the 
fluid surface ( h). These conditions are: 


w{x,y, -h B ) 
w(x,y, h,t) 


= —u 

_ dh 
~ dt 


dh B 

dx 

+ u 


dh B 

V ~dy 

dh dh 
dx dy 


(11.18) 


Combining (11.17) and (11.18) yields (after simplification outlined in Appendix A) 


dh 

dt 


d_ 

dx 


(Hu) + 


d_ 

dy 


(Hu) = 0. 


E. LINEARIZING THE SHALLOW WATER MODEL 


(11.19) 


The shallow water model in its current form is non-linear. We have three state vari¬ 
ables, u, v, and h. Each of these are defined such that u = u(x, y, t) and v = v(x , y, t) are 
the unknown fluid velocities in the x and y directions, h(x, y, t) is the unknown fluid eleva¬ 
tion above the reference level, / is the Coriolis parameter, and g is the gravity acceleration. 
We now introduce the following shorthand for partial derivatives 

„ _d_ ^ d 2 

^CL rx ") Oab rx r\j 

oa oaob 
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We wish to find a linear version of these equations. Right now, we have 


d t u + ud x u + vdyU — fv = —g d x h 

d t v + ud x v + vdyV + fu — —g d y h (11.20) 

d t h + d x (Hu) + d y (. Hv) = 0. 

Now, suppose that the bottom topography is flat such that h B is constant and u and v can 
be described by a constant mean term and a small 0(5) deviation from that value, i.e., 

u = U + u* v = V + v* H — h B + h 

To be clear, U and V are the mean velocities and Kb is the mean water elevation. Using 

these substitutions and neglecting any 0(S 2 ) terms results in the linearized form of the 

shallow water equations (see Appendix B for details): 

d t u* + Ud x u* + VdyU* - f(V + v*) = -g d x h 

d t v* + Ud x v* + Vd y v* + f(U + u*) = -g d y h (11.21) 

d t h + Ud x h + Vdyh + h B (d x u* + d y v*) = 0 . 


F. KLEIN-GORDON EQUIVALENT TO THE SHALLOW WATER MODEL 

Using the linearized form of the shallow water equations, we can find a Klein- 
Gordon equation equivalent through a series of linear operations as outlined by Pedlosky [13]. 
We begin by defining the operator (linearized Lagrangian derivative) 

^7 = dt + Ud x + Vdy. (11.22) 
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Substituting (11.22) into (11.21) we have the modified form 


Dh 

~Dt 


Du* 

Dt 

= -g d x h 

(11.23) 

Dv* 

Dt + ^ + 

— — g d y h 

(11.24) 

+ ha ( d x u + dyV ) 

= 0. 

(11.25) 


Taking / as constant (along some latitude) and summing the partial derivative of (11.23) 
with respect to x and the partial derivative of (11.24) with respect to y, we have 


o l Du * \ ff) * 
Or ( ) - fd x v 

' Dv* 

d« 1 dF 1 + fdyU * 


Q^xxh 

Q^yyh' 


( d x u* + d y v*) + / ( d y u* - d x v*) = -gS/ 2 h 


(11.26) 


Similarly, we find the difference of the partial derivative of (11.23) with respect to y and the 
partial derivative of (11.24) with respect to x to find 


d„ 


dr 


Du* 

~Dt 

Dv* 

~Dt 


fdyV* = - gd xy h 


+ fd x u* = —gd xy h 


( d y u* - d x v*) - / (, d x u * + d y v*) = 0 


(11.27) 


We apply the operator to (11.26) and add —/ times (11.27) 

^2 ( dxU * + 9 V V *) + ( d V U * ~ d * V *) = ~ 9 J^ t ( V2/l ) 

+ / 2 ( d x u* + d y v*) - f^~ t ( d y u* - d x v*) = 0 

^ (d x u* + d y v*) + f ( d x u* + d y v*) = -g^- t (V 2 h) . (11.28) 
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From (11.25), we know that = d x u* + d y v*, which can be used in (11.28) 



which, after factoring becomes 

wt& +f2h -^ v2h )=°- 

We can rewrite this equation as 

+ fh - CqV 2 /?. = S(x, y, t ) 

where f = ghs and j^(^S(x,y,t)^ = d t S + Ud x S + Vd y S = 0. This source term we 
will assume to be zero, giving us the homogeneous form 

^ + fh - c 2 0 V 2 h = 0. 

If we expand the operator twice, we get an expanded Klein-Gordon form 

(d t + Ud x + Vd y f h - CgV 2 /i + fh = 0 (11.29) 

of the shallow water equations under constant advection U and V and dispersion evidenced 
by the f term. This equation specifies the perturbation of the wave height h above a 
reference level hn- 

G. RECOVERING THE FLUID VELOCITIES 

Suppose now that we have the solution for h(x, y, t ). In order to recover the fluid 
velocities u(x,y,t ) and v(x,y,t), we consider the modified form of our shallow water 
model shown in (II.23)-(II.25). We first apply the operator f to (11.23) and multiply (11.24) 
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by / to yield: 


DV Dv* DV _ D (dh 
~dT ~ -’lJt ~ *~Dt ~ ~ 9 Dt 
Dv* 2 dh 

s ~m +f {u + u) -~ gf ai 


(11.30) 


(11.31) 


Adding (11.30) and (11.31) yields: 


wt +f2 ) u ‘ +f2u 


( D fdh\ dh 

^ \ Dt \ dx) dy 


(11.32) 


The solution of this partial differential equation (no more difficult to solve than the equation 
for the perturbation of the wave height h) gives us the fluid velocity in the x direction. 


A similar manipulation (subtracting / times (11.23) from -D operating on (11.24)) 


yields: 


^ + , 2 V + / v = - 9 M\ ^ V 


Dt \dy 


(11.33) 


The solution of this partial differential equation gives us the fluid velocity in the y direction. 
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III. HIGH ORDER NON-REFLECTING BOUNDARIES 


The numerical solution of a wave propagation problem in a very large or unbounded 
domain provides a challenging computational difficulty - namely, solving the problem on 
a finite computational domain while maintaining the true essence of the solution. One of 
the modem techniques that has garnered a significant amount of attention in handling this 
challenge is the absorbing or non-reflecting boundary condition (NRBC) method. In using 
this method, the original infinite domain is truncated by an artificial boundary B, resulting 
in a finite computational domain Q and the residual domain I). Figure 4 illustrates the 
NRBC set-up using an infinite wave guide. Here, the artificial boundary B extends from 
the southern (T$) to the northern (Tat) boundaries of the wave-guide, thus creating the east 
(r E ) and west (IV) boundaries of Q at x = x E , x w respectively. Appropriate boundary 
conditions are prescribed on the northern (Tat) and southern (T s ) boundaries. Outside of 
the area enclosed by these boundaries is the residual infinite domain D. 

y 


X 

Figure 4: An infinite wave-guide truncated by artificial boundaries Tw and F E 

One would expect that the introduction of any boundary B, where one does not 
physically exist, to pollute the solution through the reflection induced by such an artificial 
boundary. In the last two decades, significant efforts have been extended to find stable, 
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efficient, accurate and practical means of reducing this reflection through so-called NRBCs 
[18]. 

Several high-order NRBCs have been devised to reduce spurious reflections that 
would pollute the solution. Beginning in the late 1980s, the well-known Engquist-Majda 
[19] and Bayliss-Turkel conditions [1] gave way to Collino’s [2] low derivative, auxiliary 
variable formulation for the 2D scalar wave equation. This sparked a flurry of activity in an 
effort to find quality, high-order NRBCs that were easily implementable. The sheer volume 
of literature on the topic of boundary conditions for infinite problems suggests that there 
is no “perfect” boundary condition available for general purpose use. In reality, a modeler 
must make decisions on how to balance accuracy, efficiency and ease of implementation to 
yield reasonable solutions. Extensive reviews on the topic can be found in [18, 20, 21, 17] 

A. HIGDON SCHEME 

The starting point for the family of NRBCs discussed in this dissertation is the 
condition devised by Higdon in a series of papers [22, 23, 24, 25, 26, 27], that was demon¬ 
strated in a low-order finite difference setting. While in theory, Higdon’s NRBC is con¬ 
sidered a high-order NRBC, the formulation requires evaluation of increasing high-order 
spatial and temporal derivatives as the order of the NRBC is increased. Higdon’s condi¬ 
tion (and most NRBCs for that matter) seeks to annihilate waves impinging normal to a 
boundary. To see the idea behind this condition, consider a one-dimensional wave equation 

dtth Cq 0 xx h 0 

whose solution was obtained by d’Alembert in 1747 [28], as 

h(x , t ) = F(x — c 0 t) + G(x + c 0 t). 

This solution implies that there are two components to the solution - one wave G of fixed 
shape moving to the left at velocity —cq and one wave F of fixed shape moving to the 


24 



right at velocity c 0 . Now, suppose that the right moving wave approaches a boundary. To 
perfectly absorb the wave impinging on the boundary, the boundary must satisfy G = 0 
such that the boundary condition is h(x,t ) = F(x — Cat). Differentiating the boundary 
condition with respect to x and t results in: 

d x h = F'{x — cot), d t h = — c 0 F\x — c 0 t), (III. 1) 


which implies 


d t h + c 0 d x h = 0. 


(HI-2) 


This is the Sommerfeld radiation condition for the eastern boundary. If we expand the 
discussion to two-dimensional problems, this condition implicitly assumes that by the time 
the wave front reaches the boundary, it is traveling primarily as a plane wave at speed c 0 . 


1. Accounting for Dispersion 

In a non-dispersive medium, waves can propagate without deformation. The chal¬ 
lenge associated with dispersive waves such as the Klein-Gordon equivalent under consid¬ 
eration here, is that the wave speed depends on the frequency of the wave. Thus, if using the 
Sommerfeld radiation condition (III.2), only the waves traveling at phase speed c 0 will be 
absorbed - for all others, only a portion of the wave will be absorbed. Higdon considered 
a composition of J of these simple first order operators to yield his boundary condition: 


Hj: 




h = 0 


on T 


(HI-3) 


where d n is the normal derivative on the boundary. In the case of the wave guide shown in 
Figure 4, this derivative is d x and —d x corresponding to the eastern and western boundaries 
respectively. The boundary condition contains parameters Cj that can be interpreted in 
terms of the phase velocities of waves absorbed exactly at the boundary. Except in contrived 
examples, there are infinitely many waves composing the solution, and in a dispersive 
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medium, a corresponding infinitely many phase velocities. A choice of the order J of the 
boundary condition seeks to annihilate the “most significant” J waves. 

2. Reflection Caused by the Boundary 

For purposes of illustration, consider a simplified version of (11.29) where U, V are 
both zero 

d tt h - c 2 0 V 2 h + f 2 h = 0. (III.4) 

Further, suppose that the domain is structured such that the NRBC is imposed on only on 
the eastern boundary of a rectangular grid as shown in Figure 5. On the south and north 


r N 



Figure 5: A semi-infinite channel truncated by artificial boundary Te 
boundaries (r^ and T N ), we have no normal flow, i.e.,: 

d y h = 0 on Tat and (III.5) 

and we impose h = f(y, t) on Tw- As x —> oo the solution is known to be bounded and 
not to include any incoming waves. A solution to (III.4) has the form 

h(x, y, t ) = Y(y) cos ( kx — c ot + 4>) (III.6) 

such that Y(y) satisfies (III.4). One such example Y(y) = cos ^ for n — 0,1,2,... that 
satisfies these boundary conditions is given by Givoli and Neta [29]. Given this choice for 
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Y(y), the dispersion relation is 


w 2 = cl (V + + f, n — 0,1,2,.... (IIL7) 

In this solution, k is the horizontal wavenumber, n is a parameter for controlling the shape 
of Y(y), lo is the wave frequency and / is the phase shift. The horizontal phase velocity [28] 
is therefore C k = f for a particular wave number. Suppose that one of the C/s in (III.3) 
equals C k . 


d t h = uY(y ) sin (kx — ut + /) 

d x h = —kY'{y) sin {kx — c ut + 4>) 
to 

so that 0 = d t h + — d x h 

k 

thus, satisfying the Higdon boundary condition (III.3) exactly for that particular mode. 
If, however, none of the C/s were identically C k , then a portion of the mode would be 
reflected and the boundary condition would not be exactly satisfied. To make the boundary 
condition true, it would have to be adjusted to incorporate the reflected modes. Higdon [27] 
sketches and Dea [17] details via a simple induction argument that proves the reflection 
coefficient Rj is 

*=ri 

j =i 

We notice here that Rj is a product of factors that are each less than one. Therefore, 
simply increasing the order of the NRBC (J) reduces the amplitude of the reflected wave 
irrespective of the choice of C r van Joolen et al. [30, page 1045] notes, “Of course, a 
good choice for the C 3 would lead to better accuracy with a lower order J, but even if the 
‘wrong’ 6/s are taken .. .one is still guaranteed to reduce the spurious reflection as the 
order J increases.” We can use the dispersion relation together with the definition of the 


C 3 - C k 


Ci + c k 


(HI-8) 
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horizontal phase velocity C k to find that 


C k 


( k + W )+/ 2 

k 



2 nV , f 2 
L 0 ~r~ J 

k 2 


which shows that C k > c 0 and therefore guides our selection of C 3 to always be at least c 0 . 

Based on the preceding discussion, we outline some of the inviting characteristics 
of the Higdon NRBC. 

• The boundary condition is tractable, extending basic principles to arrive at the final 
condition. 

• They are exact for all waves that propagate with speed C r 

• Reflection is guaranteed to decrease by simply increasing the order of the NRBC. 

• They have been applied to a wide variety of wave-type problems including those in 
dispersive media [17, 22, 23, 24, 25, 26, 27, 29, 31, 32]. 

The boundary condition, however, suffers from an implementation point of view. Due to 
the increasing order in the spatial and temporal derivatives required, even Higdon in his 
original papers considered practical implementation to be no more than J = 3. 


B. HIGDON ADJUSTMENTS FOR ADVECTION 


In the discussion of Higdon’s boundary condition above, we considered the zero 
advection case. Eventually we wish to implement the Klein-Gordon equivalent which in¬ 
cludes constant advection. Here, we discuss modifications to Higdon’s scheme to accom¬ 
modate non-zero, constant advection. This discussion will use physical arguments that will 
be reinforced by numerical considerations later in this dissertation. To begin, we consider a 
wave augmented by advection impinging on the NRBC of the semi-infinite channel shown 
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in Figure 5. Implicit in this derivation is that by the time the wave pulse arrives at the 
NRBC, it is traveling primarily as a plane wave (in the x— direction). The wave moves 
according to the equation under consideration, 

(d t + Ud x + Vd y f h - c 2 0 V 2 h + f 2 h = 0. 

If, however, this wave is moving primarily in the a;—direction, then any effects of y are 
negligible. Further, as in the derivation of the Sommerfeld radiation condition (III.2), we 
first consider a non-dispersive environment such that f 2 = 0 in this derivation to suggest 
an appropriate boundary condition. This simplifies our discussion to a wave that moves 
according to 

(d t + Ud x f h - c 2 0 d xx h = 0. (HI.9) 

As outlined in Appendix C, the solution takes the form 

h(x,t ) = F^x — (c 0 + U)t^j + G^x + (co — U)tj (III. 10) 

with the interpretation that the general solution is the sum of F, a wave of fixed shape 
moving to the right with velocity c 0 + U and G, a wave of fixed shape moving to the left 
with velocity c 0 — U. 

As described in III.A, if we consider the wave moving to the right approaching the 
eastern boundary, the boundary must satisfy G — 0 such that the boundary condition is 
h(x, t ) = F (x — (c 0 + Ujt'j. Differentiating the boundary condition with respect to x and 
t results in: 

d x h = F’[x- (c 0 + U)t) dth = -(co + U)F’(x - (co + U)tj, 

which implies 

d t h + (cq + U)d x h = 0. 
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Following the discussion for the reflection coefficient, this suggests the “educated” choice 
of augmenting the Higdon parameter Cj with the added advection. This choice will be 
revisited later in this dissertation. 

C. GIVOLI-NETA AUXILIARY VARIABLE LORMULATION 


In [29], Givoli and Neta directly extended the Higdon scheme to high-order finite 
difference discretizations via an algorithm where the order of the NRBC was simply an 
input parameter. They later extended this formulation [33] to one that does not involve 
any high derivatives (hereafter referred to as the G-N formulation). The elimination of all 
high-order derivatives is enabled through the introduction of special auxiliary variables on 
B. This construction demonstrated in [33] and [34] for finite differences was further ex¬ 
tended in [35] for finite element schemes to solve the dispersive wave equation. Hagstrom 
and Warburton [36] also used the Higdon and auxiliary variable framework to develop a 
symmetric boundary formulation in a full-space configuration where special comer com¬ 
patibility conditions were developed for the non-dispersive wave equation. Extensions and 
comparisons between the two methods were published by Givoli and Hagstrom et al. in 
[37] and [38], 

We present a brief summary of the G-N auxiliary variable process as described in 
[35]. For the semi-infinite channel shown in Figure 5, this auxiliary formulation begins 
with the Higdon boundary condition given by: 


Hj : 



(III.ll) 


Auxiliary functions <j)\, ..., 4>j~ i, which are defined on Tg as well as in the exterior do¬ 
main D are now introduced. Eventually, we shall use these functions only on Ye, but the 
derivation requires that they be defined in D as well, or at least in a non-vanishing region 
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adjacent to T e- The functions <j)j are defined via the relations 


{ a * + k a ‘) 


h = 0 i , 


(III. 12) 
(III. 13) 


(d x +-^dt^ <l>j-i = 0. (III. 14) 

By definition, these relations hold in D, and also on Ye- It is easy to see that (III.12 - 
III. 14), when imposed as boundary conditions on Ye, are equivalent to the single boundary 
condition (III. 11). If we also define 


0o = h 0j = 0 , 


(III. 15) 


then we can write (III.12 - III.14) concisely as 

(d x + ^<9^ 0j_i = 0J j = (III. 16) 

This set of conditions involves only first-order derivatives. However, due to the appearance 

of the ^-derivative in (III. 16), one cannot discretize the <i> 3 on the boundary T E alone. 

Therefore we shall manipulate (III. 16) in order to get rid of the ^-derivative. 

The function h satisfies the dispersive, advective wave equation (11.29) in I). Since 

the function 0! is obtained by applying the linear (constant coefficient) operator {d x + ~^rdt 

to h, it is can be shown that 0 \ should also satisfy the same equation in I) 4 . Further, since 

cj)j is obtained by applying the same linear operator j — 1 times to 0 1 , the functions 0, 

4 Here we must use the assumption that c f) and f are constants. By applying the differential operator to 
(11.29), computing each of the <j)j derivatives present in (III. 17) using the differential operator and simplifying, 
a simple induction argument shows that the cf>j 's must satisfy (III. 17) 
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should satisfy an equation like (11.29), namely, 


(dtt + (U~ ~ Cq) d xx + {V~ — eg) 7 y . y + 

2(7 <9^ + 2(7 d y t + 2(717 <9 xy + / 2 ^ 77 = 0 (III. 17) 

Using the following identities: 

d X x4>j = (&x ~ <9* j ( 9 X + T7 dt j (f)j + —^2 4>j 

V °j+i / V W+i 7 L y+i 

dxt&j 7/ (7,x 0j ) 

7;/:y77 7y (7 x4>j) 


and combining with (III.16) allows us to write (III.17) as: 


/2f7 (7 2 -c 2 \ 

U C 2 J 


4-i + 


2(717 



- 21") - (l" 2 - c 2 ) 

- 2t/) 7 - 2 UV# S - Z 2 ^-1 = 


(U 2 - e 2 ) 0 J+1 


for j = (III. 18) 


The details of this transformation are shown in Appendix D. In (III. 18) and elsewhere, a 
prime indicates differentiation with respect to y along T E , be., the tangential derivative 
along T E . As desired, the new boundary condition (III. 18) does not involve ^-derivatives. 
In addition, there are no high-?/ or t derivatives beyond second order. It should be noted that 
van Joolen et al. [31] developed an equivalent formulation in using Lagrangian derivatives. 
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Rewriting (III. 12), (III. 18) and (III. 15), the new formulation of the J^-order NRBC 
onTg can be summarized as follows: 


/3 0 h + d x h = , 

a j4>j -1 + ~~ \ftj -1 + Pj&j ~ 7 ' f 2( Pj -1 = 


= h 


0j = 0 


(III. 19) 
(III. 20) 
(III.21) 


where 


A) 

& 


5 q 


1 21/ U 2 — ci 

- o • =-1-- 

r, ’ ■' Cj c) ’ 

(U 2 - c 2 ) (7 + 7-) _ 2(7. 7 = 2W, 


«< = 77 _ 2V, A, = l/ 2 - cS, 


A 7 . = U 2 -cl 


What remains is to link the boundary condition to the interior formulation. As will be 
shown in the next chapter, the means to link the two formulations will naturally follow 
from the spectral element framework. 
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IV. DISCRETIZATION VIA SPECTRAL ELEMENTS 


Once a suitable NRBC is devised, the problem must be solved numerically in O by 
the finite difference, or, as in the case of this analysis, the spectral element (SE) method. 
The SE method, originally introduced by Patera, .. combines the generality of the finite 
element method with the accuracy of spectral techniques ...” [39, p. 468]. The key to the 
spectral element (SE) method is the careful selection of the integration and interpolation 
points in order to yield accurate but efficient solutions. 

As indicated in Chapter III, the Givoli-Neta auxiliary variable formulation has been 
previously demonstrated in both finite difference and finite element schemes to arbitrarily 
high NRBC order, however, accuracy gains realized by increasing the NRBC order slowed 
significantly after certain orders. The SE formulation used in this dissertation seeks to rem¬ 
edy this limitation by using a high-order treatment of space (SE) to show the benefits of 
using a high-order boundary (G-N) scheme. It should be noted that the only other spec¬ 
tral element, high-order boundary approach (to the author’s knowledge) is [40] that is in 
press. That paper shows how a SE approach combined with high-order boundary treatment 
significantly improves the accuracy for the non-dispersive wave equation on a semi-infinite 
channel using the Hagstrom-Warburton formulation [36]. The key difference in our work is 
that we use high order space, boundary and time integration (discussed later in this disser¬ 
tation) in both a dispersive and non-dispersive wave equation setting. We show the details 
of early results (absent of advection) in [41]. 

What follows in this chapter is a brief overview of the SE method as presented by 
Giraldo [42]. Additional treatment of this method can be found in Giraldo-Restelli [43] or 
Pozrikidis [44]. For this problem, we will discuss two formulations - one for the interior 
and one for implementing the boundary conditions. 
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A. 


INTERIOR WEAK INTEGRAL FORMULATION 


Finite difference schemes seek to discretize differential operators based on Taylor 
series representations of the function near a point. Spectral element schemes seek to solve 
the equations in a weak integral form. It can be shown that if the governing equations 
satisfy the weak integral form, then they also satisfy the differential form provided the test 
functions used are sufficiently differentiable. 

The weak form of the problem in Q is now constructed for the semi-infinite channel 
described in III.A.2. Following a standard Galerkin approach, the solution is sought in the 
space of test functions, 

V = {h\h E H\ ft) and h = f(y,t) on T w }. (IV. 1) 


Here, is the Sobolev space of functions whose first derivatives are also 

square-integrable in fL Now, Equation (11.29) is multiplied by the globally defined basis 
functions T y) G V and integrated over O so the weak form is: 


r .. r d 2 h 

/ * i hdn + X x / T, 

In Jn c,x 


(IQ + Xy / ^ 


d 2 h 
dy 2 


dQ 


+ 2 U I ^~dn + 2V [ 

Jn (xx Jq dy 

+UV I T, ( t) ' h ■ n ' h ) dQ + f 2 [ *ih dn. = 0. 
Jn \oxoy oyoxj Jq 


(IV.2) 


Here, h and h are the first second derivatives of h with respect to time. Further, || and 
are the mixed second derivatives of h with respect to t , x and y. 

In order to maintain a symmetric form of the problem, the mixed derivative has 
been split appropriately. To ensure the solution h is in the vector space // 1 requires some 
special handling of the second order derivatives, which is facilitated by the use of Green’s 
theorem, i.e., 


T Q2k .70 

^. t7 — dU = 
ox- 


T dh ,-r^ 
dr - 

ox 


d^i dh 
dx dx 


dn 
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where n is the outward normal on Y and n x is the ^-component of that outward normal. 


Extending this idea for each second order derivative in (IV.2) gives the weak form: 


/‘ t v x /' d^i dh ^ r d^i dh _ 

/ \l hh dQ- X x / -7— 7T- dU — X y / -7—77- dtt 

Jn Jn dx dx J n dy dy 

-uv l + dn 

Jn \ dy ox ox dy J 

+ 2 U I v~dn + 2 V I v—dn + f I r t hdn (iv.3) 

Jn dx Jn ay Jn 

f dh f dh f (dh dh \ 

+Aa; / ^i—n x dr + \y / dT + UV / 'k l —n y + — n x dr = 0. 

Jr ox y J r dy J v \dx dy J 


Recall that for the semi-infinite channel described in Figure 5, the corresponding 
boundary conditions for T N andT s are d y h = 0, displacement is specified on and 
Ts is a non-reflecting boundary. Using these boundary conditions on the northern and 
southern borders along with the normal vectors specified by the structured, rectangular 
grid shown in Figure 6, the problem may be simplified to account for contributions on 
individual boundaries. Using this information along with (III. 19) to eliminate the normal 
derivative term on Tg, we get: 



Figure 6: Normal Derivatives to Boundaries 
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in 


T V 9h . 

'Ifi/i o?n — A rr / ——— o?n — A,, 

or ox 


<9\l/, <9/r 


i/x i/x dy dy 

_ uy f fd%dh + d%dh 

Jn V d y ® x dx ®y 


dQ 

dn 


f f)h f f)h C 

+2U / dn + 2v / Vi— dn + f / %h<m 

Jn dx Jn dy Jn 


+A a 


The)i cir i 


do 


Tq/i cir e 


/> ai /» o 7 

+uv / rfr^v -uv / ^7T rfr s + uv. *,» 

arv ar s ^ Jr E dy 


r E 

T <9/r 

\n7~ cir£ = 0. 


(IV.4) 


Note, that two of the boundary integrals (other than those on T e) survive after simplification 
- namely on T^ and Y s ■ This occurs since the boundary condition (III.5) only specifies no 
flux (i.e., reflecting boundary conditions) on these boundaries. This is a common boundary 
condition for inviscid flow problems. 5 


B. BOUNDARY WEAK INTEGRAL FORMULATION 


Since the term <J>i appears in the interior formulation (IV.4), this is not a complete 
weak form of the problem in Q. We turn our attention to computing a separate weak form 
for (III.20) to complete the problem statement. This solution is sought in the space of test 
functions, 

Vr E = {</>j\<t> j eH 1 (T E )} (IV.5) 

where II 1 (r E ) is the Sobelov space of functions whose first derivatives are also square 
integrable onT B . 

As in the interior formulation, we multiply (III.20) by the globally defined, one¬ 
dimensional basis functions (on the boundary) Q e Vr E and integrate it over T e- After 

5 It should be noted that in this analysis, the two-dimensional basis functions, 'It,, are constructed using 
tensor products of one-dimensional basis functions, Q. This means that on the boundaries, these T, basis 
functions collapse to Q on the boundary. Therefore, in practice the boundary integrals are constructed using 
Q basis functions, requiring far less storage. 
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integration by parts and simplifying this yields: 


a,- 


Ci4>j-1 d^E + Kj / 007 — 1 dT E + \ y / Ciftj- 1 + f3j / 


- 7 / (iftj dr E - f 2 dr E = \ x / C:7j+i 


(iv.6) 


for j = 1,..., J — 1. As in the interior formulation (IV.2), 0 and 0 are the first and second 
derivatives of <f> with respect to time and qY is the tangential time derivative of (f). Recall 
from the auxiliary variable formulation (III.21) that 0 O = h and (f)j = 0. 

The formal problem statement is then: Find h e V and 0 3 e W> ; where j = 
1,..., J — 1, such that Equations (IV.4) and (IV.6) are satisfied V T, e V and Q E VV ;; . 

C. GRID GENERATION AND CHOICE OF BASIS FUNCTIONS 


One of the key advantages of the spectral (finite) element formulation over differ¬ 
encing schemes is the ability to represent complex geometries by breaking up the domain 
into simple geometric shapes. Triangles and quadrilaterals are primarily used to represent 
these geometries in two-dimensions, while tetrahedra and hexahedra are primarily used to 
represent geometries in three-dimensions. The sheer volume of grid generation software 
available for generating triangles (tetrahedra) and the dearth of grid generation software for 
quadrilaterals (hexahedra) suggests that the former are the more common means of repre¬ 
senting the geometry of a problem. However, there are some very nice properties of the 
latter that led to the choice of structured (and later unstructured) quadrilaterals to discretize 
the geometry in this analysis. 

While triangular grids are typically easier to generate, quadrilateral elements can 
be formed directly, or more commonly, by combining two or four basic triangle elements 
as shown in Figure 7 6 . This idea can then be used to find an appropriate mesh for simple 
or complex geometries. For this analysis, unstructured quadrilaterals were generated using 
6 Adapted from [45], Figure 5.2, p. 141. 
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Figure 7: Arbitrary quadrilateral element formed by combining triangles. 


Automesh-2D [46], an automatic mesh generator that provides high-quality results for the 
domains considered in this dissertation. A sample of the mesh that can be generated is 
shown in Figure 8 7 . Noteworthy in this mesh is that the complex geometry of each letter 
has been broken up into approximately 200 non-overlapping quadrilateral elements Q e 
such that 

Nett 200 

o = (J o e . 

e=l 

Further, each Vt e is constructed such that all internal angles are all between 30° and 150° - 
a key criterion in providing a stable numerical solution. 



Figure 8: Sample mesh generated using Automesh-2D. 

To see why this criterion is important, consider a mapping from the physical co¬ 
ordinate space to a canonical element space such that x = x(£, //) and y = y(£, rf) where 
7 Geometry for the letters graciously provided by [46]. 







(£, rj) e tt c = [—1, l] 2 . This nonsingular mapping is made to facilitate efficient and accu¬ 
rate computation of operations such as differentiation and integration [47]. One such ele¬ 
ment that illustrates this mapping can be seen in Figure 9. All derivatives are then mapped 


y 




Figure 9: Mapping from physical space to computational space 


to the canonical element space by virtue of the chain rule to yield the total derivative: 

( dx\ / if ip \ / df \ 

UHtsJU)' 

From the total derivative, we can find key components required for the evaluation 
of the integrals induced by this canonical element mapping [42]. Such terms include the 
transformation Jacobian and its associated determinant: 


( dx dx \ 

an a v 

Oy Oy J ’ 

dt dy J 


dx dy dy dx 
d£ dy dC, dy 1 


as well as mappings of metric terms: 


(IV.7) 


d £ 1 dy <9£ 1 dx dy 

dx | J e | dy' dy \J e \ dy ’ dx 


1 dy dy 1 dx 
J e \d £’ dy | J e | <9£ ’ 


(IV.8) 


If any of the nodes comprising the elements approach interior angles of 180° (degenerating 
the quadrilateral into a triangle), it can be shown that the transformation Jacobian at that 
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point tends to 0 causing each of the metric terms to tend toward infinity thus corrupting the 
solution. Derivation of the metric term mappings and consequences of degeneration of the 
quadrilateral are detailed in Appendix E. 

Crucial to the construction of SE approximations is the representation of the solu¬ 
tion variables in terms of smooth basis functions 'Iq. While in theory these basis functions 
are defined throughout the entire domain they have compact support, and are in practice 
constructed locally. This implies that we can solve the global problem by simply summing 
up the contributions from the smaller elemental problems in a process known as direct 
stiffness summation. Ideally, these basis functions will not only support the grid chosen, 
but will also have properties that facilitate the numerical integration of the weak formula¬ 
tion [48]. 

For this analysis, the local element-wise solution h is represented by N lh order 
Lagrange polynomials in (£, rj) such that 

M N 

v) = (£’ *7) h *7*) 

k= 1 

where (ff, rj k ) are the M N = (N + l) 2 Legendre-Gauss-Lobatto (LGL) interpolation points 
and iff are the associated multivariate Lagrange polynomial basis functions. The square 
structure of the canonical element simplifies matters in that we can express iff by a tensor- 
product of the one-dimensional Lagrange polynomial basis functions as 

k (€, v) = u i (0 Vj iff) 

where iff = 1, 2,.... A r + 1 and k = 1, 2,..., M N . Further, //, and i/j are the one¬ 
dimensional basis functions generated using the LGL points in £ and //. To get from the 
one-dimensional local indices (iff) to the two-dimensional local index k requires the map¬ 
ping k = i + (j -1)(N + 1). 
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The choice of the LGL points for interpolant generation have the especially attrac¬ 
tive property of allowing efficient quadrature. In particular, any polynomial of up to degree 
2 N — 1 can be integrated exactly (to machine precision) by evaluating it at only N + 1 
points. Additionally, if the interpolation points are also used as the sampling points in the 
Gaussian quadrature rule, it has the added benefit of yielding a diagonal mass matrix as the 
LGL-based basis functions 'F satisfy a discrete orthogonality via the cardinality property 
of the Lagrange basis functions. The effect of these choices is to allow the appropriately 
weighted summation of the basis function expansion of h to evaluate all integrals. There¬ 
fore, to approximate a local elemental integral, we find 



h(x, y)dfl e 


h(t,-n)\j e \dn c 


N +1 

~ Y u (&) u fai) h (&> Vj) \ je (&> Vj) 

i,j= 1 


where u (£*) and u (■ r/j ) are the quadrature weights associated with one-dimensional LGL 
quadrature and is the canonical element region of integration. This process is described 
in greater detail in [42]. 


D. GALERKIN EXPANSION 


We now turn our attention to the spatial discretization. First, we construct the N th - 
order approximation of the variables h and (j>j using the same basis functions used in the 
weak form IV.4 and IV.6 as follows: 

N P N b 

h N = J2 *k hk > = Y, J = 1,2,..., J-l. 

k= 1 k =1 

Here, N p and A), are the number of points that O and T e are discretized into, repsectively. 
For a structured quadrilateral grid with A(f and AT| elements in the x and y directions, 
N p = (N^N + 1) (N%N + 1) and Nb = N^N + 1. Next, we substitute this basis function 
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expansion directly into the weak formulations, resulting in the following interior 


Mh + (2 U D x + 2V Dy) h + (—A X L XX — A y L yy + f 2 M — UV (L yx + L xy )) h (IV.9) 
+\ x A E <t) l - p 0 \ x B E h + (UVC N - UVC S + UVC E ) h = 0 


and boundary 


ajM %_i + CjD%-! + (A y L b - f 2 M b ) (IV. 10) 

+l3 J M b d>j - 7 D b (f>j - \ x M b (j) j+ i = 0, j = 1,..., J - 1 


forms of the problem. Here, M, D x , D y , L xx , L yy , L xy and L yx are interior formulation 
matrices of size N p x N p . Further, B E , C E , C$ and C E are interior formulation matrices 
integrated along the boundaries of size N p x N p while A E is of size N p x N b . Finally, 
M b , I) 1 ' and L b are auxiliary formulation matrices of size N b x N b . These global matrices 
are obtained from analogous element arrays via assembly, given by: 


N e 


N e 


M=/\ M e , B E =/\ B%, 

e=l e=l 

N e N e 

L X X = /\ L xx , Lyy = L yy , 

e=l e=l 

N b N b 

M b = /\ m b , D b = /\ d b , 

6=1 6=1 

N e N e 

c s = /\ c%, c E = /\ c%, 


e=l 


e=l 


N e 


N e 


d* = A a > = A D l 

e=l e=l 

N e N e 

Lxy = f\ L xy , L yx = L yx , 

e=l e=l 

N b N e 

l" = A i\ c« = A °N’ 

6=1 

N e 

Ae — f \ A e E 


e=l 


e=l 
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N e 

where j\ is the direct stiffness summation operator required by all continuous Galerkin 

e=l 

methods. The expressions for the element arrays are: 


M e - = 
13 


D e ■ = 
y,v 


r e = 

xy,ij 


dt = 


s~ie _ 

U s,ij - 


i’ii’j dfle 




f dipt dtfjj 


dy dx 


v r v' s dF e 


dQ e 
dQ P 


1 r§ 


$ 


rf r s 

■ dx ■ 


ne _ 

n E,ij ~ 


T e — 
xx,ij 


T e = 

yx,ij 


l = 

v rr* C 


s~ie _ 

U E,ij - 


I ipiipj dT e 

/ O (Xu Ug 

f£l e dx OX 

d ^ 7 T da ' 

In. dx ay 


v'Xs dF e 


IpE 




di± dT E 

dy dl ‘ 


D e = 

X,IJ 


T e = 
yy,ij 


m rs = 


r u _ 

^N,ij — 


A e 

EAr ~ 


f ipzjp- dQ 

I n e dx 

f d'lpi dipj 


In e dy dy 
[ 13 r i3 s dr e 


dfL 


IpN 


ibAl dr N 

V'dx ‘ 


I P E 


ipii3 r dr 


where fl e and T e denote, the part of Q and T associated with element e. Also, VJ t and u, 
are the locally defined basis functions from which the global basis functions (4q and Q) are 
formed. For quadrilateral elements with polynomial order N as used in this analysis, ip % will 
be discretized into (N +1) 2 points, and u t into (N +1) points, thus, i,j = 1,2,... (A^ +1) 2 
and r, s — 1, 2,... N + 1. 

Now, let: 


A = M" 1 (A,A x B e - 2 UD X - 2 VD y ) 

1 = M - 1 (A X L X + XyL y - f 2 M + UV (L xy + L yx - C N + C s - C E )) (IV. 11) 

C = -X x M~ 1 A b E 


and substitute equations (IV. 11) into (IV.9) to get the matrix form of the interior problem: 


h = Ah + Mh + C(f>i 


(IV. 12) 
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If we examine the boundary auxiliary variable formulation (III.20), we see that 
the selection of appropriate Cj values for the auxiliary variables has not yet been fully 
addressed. As discussed in III.A.2, any choice of Cj is guaranteed to reduce spurious 
reflection as the order of the NRBC (J) increases, however based on the discussion in 
III.B, the Higdon parameter should somehow be augmented with the advection. 

Armed with this, we choose convenient values for our C/s that cause the second 
order in time (ay) terms to vanish. In the case of the semi-infinite channel, on the eastern 
boundary this value is: Cj — c 0 + U. A physical argument for this choice is that since the 
predominant speed of the wave, absent the advection terms, is cq, the Cj terms “correct” the 
boundary formulation to account for the advection. This selection for the C'/s then makes: 

oli — oli — --- — CX-J -1 = 0 , 

«1 = K 2 = ... = Kj -1 = K, 

Pi — P 2 = ■ ■ ■ — Pj -1 = P 

Now, let: 

D = f 2 M b — \ y L b (IV. 13) 

and substitute (IV. 13) into (IV. 10) to get the following form of the problem: 

PM b j >i = D h rE -KD b ii r + 7 D b fa + \ x M b fa 
k D b fa + pM b fa = D fa + 7 D b fa + A x M b fa 
k D b fa + /3M b fa = Bfa + yD b fa + X x M b fa 

K D b fa-2 + PM b fa-i = D fa-2 + 7 D b fa _! 
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where h-p E is the value of h on If we now collect the terms on the left and right, we get 
the matrix form of the problem: 

E<f> = F<f> + /i (IV. 14) 

where: 
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V. SOLUTION OF THE DYNAMIC SYSTEM 


The equations described by (IV. 12) and (IV. 14) constitute a coupled system of 
ODEs that must be solved for h(x,y,t). Our approach uses standard k th order explicit 
Runge-Kutta (RK) methods to integrate the system in time. Recall that RK is used to solve 
first order ODEs, and as such, the second order systems described must be converted into 
a larger system of first order ODEs. Using the substitution v — h , v — h yields the first 
order systems: 


/ 0 


h 


0 / 


h 


0 




— 




+ 


0 / 


V 


1 A 


V 


C01 


E<i> = F4> + h. 


(V.l) 


This system was solved using a two-stage approach at each time step. First, the auxiliary 
system was solved to find the component <pi required for the main system. Then the main 
system was solved to find h at the next time step. This process was continued until the final 
time t f was reached. 

A. RUNGE-KUTTA METHODS 


RK methods are convenient in this analysis in that the machinery to decrease the 
truncation error is the same for 2 nd , 3 rd or even k th order schemes. While there are other 
schemes that require fewer function evaluations and, in fact, are more efficient than the 
explicit RK schemes used in this dissertation, we desired a method that could change the 
order of the time integration method simply as a parameter while using the same basic 
formulation. RK methods use a standard formulation that is outlined as follows. 
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To obtain an approximation to the solution of the initial value problem 


y' = f (t,y(t)), y{t 0 ) = yo 

a successive approximation y n+1 to y n is given by 

S 

y n +1 = y n + At y, biki 


1=1 


where k % — f ( t n + qAf, y n + At ^ a^kj ) i — 1,2,... s. 

V j=i 

Here, the coefficients , b t and c, are given by the Butcher tableau for the given RK 
scheme. The classical explicit RK-4 scheme is given by 


Cl 

On 

«12 

a 13 

a 44 

0 

0 

0 

0 

0 

C 2 

«21 

C'22 

®23 

«24 

1 

2 

1 

2 

0 

0 

0 

C 3 

«31 

«32 

«33 

034 

=*• 1 

0 

1 

2 

0 

0 

c 4 

0 4 i 

0 4 2 

a 43 

Q 44 

1 

0 

0 

1 

0 


bi 

&2 

h 

64 


1 

6 

1 

3 

1 

3 

1 

6 


For this analysis, all coefficients were computed using Maple™ using routines coded by 
Stone [49]. These routines allow for high precision computation of the coefficients re¬ 
quired for each high-order tableau. Of course, the computational complexity increases for 
increased RK order. To be specific, the required number of function evaluations (stages) 
and non-zero parameters required for the k th order RK method are given in Table 1 

Table 1: Computational cost for RK methods used. 


RK Order 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Stages 

2 

3 

4 

6 

7 

9 

11 

16 

18 

Parameters 

3 

8 

10 

22 

33 

45 

59 

102 

113 
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B. SIMPLE EXAMPLE OF CONVERGENCE 


To see the effect of high-order time integration, consider the simple initial value 
problem 

^ = f(t,y) = (2-t)y,. y{ 0) = e“ 2 (V.2) 

whose solution is 

y(t) = e _0 - 5 ( t_2 ) 2 . 


Now, we integrate (V.2) in time using our RK schemes (RK-2 - RK-10) from t — 0 to t — 5 
and quantify the errors observed over time between the exact and RK solution using the 
normalized L 2 error defined as follows: 


error\ \ i? 


( 


N 


n= 1 


^2 (Vn ~ y(nAt +1 0 ) 


2 \ 


1 

2 


^y(nAt + t 0 )' 2 

n =1 


where y n and y(nAt + t 0 ) are the numerical and reference solutions at a given time. 

Using a time-step At = 0.05, we see exponential convergence (to near machine 
precision) of the error as indicated in Figure 10. It should be noted that in this example, 
even if an extremely fine time step is chosen, the computed error using lower order RK 
methods cannot achieve machine precision. 

Take for instance the RK-4 method that has a truncation error of O (At 4 ). A step 
size of At = 10 -4 should result in an error O (10 -16 ), yet experimentally, the error is 
O (10 -12 ). This apparent discrepancy can be explained by the round-off error due to the 
finite-precision arithmetic (double) used in the computation. This error increases with the 
total number of integration steps used. This implies that to integrate (V.2) from t — 0 
to t = 5, we must evaluate f(t,y) a total of 200,000 times for RK-4 at a time-step of 
At = 10 -4 . If we compare this with a time-step of At = 0.025 using RK-10, this should 
result in an error O (10 -16 ). To integrate (V.2) using RK-10 for the same period of time, 
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Figure 10: L 2 error in solution of (V.2) using RK-2 - RK-10 time integration. 


we must evaluate f(t, y ) only 3,600 times. Experimentally, RK-10 achieves this O (10 -16 ) 
error convergence, thus highlighting the importance of round-off error in the convergence 
of the time integration scheme. 

To check the convergence of these RK methods for reasonable time-step values 
At G [0.001,1], consider the rate of convergence adapted from [50] and defined as follows: 


rate = 


log (error Atn+1 ) - log (error Atn ) 
log (At n+ i) - log (At n ) 


(V.3) 


This rate is a measure of how rapidly a given time integration method converges as a func¬ 
tion of the time-step refinement. Figure 11 shows the normalized L 2 error as a function 
of the time-step refinement. For each RK order, the rate of convergence is averaged over 
all At refinement levels. This average rate is annotated for each RK order and shows, as 
expected, that the maximum rate of convergence is the RK order. From this figure, we 
see that this theoretical rate of convergence is nearly reached in every case. It should be 
noted that once errors reach O (10 -15 ), round off errors prevent further improvements us¬ 
ing a given RK method, and as already discussed, these errors actually prevent theoretical 
improvements of lower order (RK-2 and RK-4) methods. 
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Figure 11: L 2 error as a function of time-step refinement. Rates as defined by (V.3) should 
correspond to RK order. 


Of course, one might wonder why such high precision when integrating in time 
would ever be required. For most engineering applications, where the model parameters 
themselves are approximated, this argument is a strong one. Perhaps a lower order time 
integrator would result in errors that are on the same order as the model parameters and 
any computational effort incurred in using a high-order time integrator would be wasted. 
However, since this work centers on the use of high-order spectral elements along with 
high-order boundary conditions; high-order time-integrators were also explored to examine 
all of the limiting factors of high-order accuracy. 

We propose that in order to see the full improvements of high-order boundary con¬ 
ditions requires a balance of truncation errors (and round-off errors) between all of the 
components of the numerical model; this includes the boundary conditions, the spatial 
discretization method, and the time-integrators. Practically speaking, we believe that the 
mathematical formulation should be the strong point of any model evaluation, specifically, 
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the order of the method should be chosen to match (or exceed) other errors inherent in the 


model. Experiments in this dissertation were conducted using boundary conditions up to 
order 20, SE polynomials up to order 16, and time-integrators up to order 10. 
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VI. NUMERICAL EXPERIMENTS 


Several numerical experiments were performed to solve the Klein-Gordon equiv¬ 
alent (11.29) with and without dispersion and advection. Additionally, we consider five 
primary configurations: semi-infinite channel, infinite channel, quarter plane, open plane, 
and open unstructured plane. Each of these configurations will be described in detail later 
in this chapter. 

In order to simplify the numerical simulation of the problem at hand, the KGE 
(11.29) is converted to anon-dimensional form as described in [51]. Using typical mesoscales 
in the ocean, the length scales were chosen 0(100 km), vertical depth scales 0(100 m) and 
the dispersion parameter / for Coriolis O(10 -4 s -1 ). Majda [51, page 61] describes typical 
advective velocities as roughly yjk that of the medium wave speed c 0 . Here, we choose 
to challenge the boundary further by choosing faster advective velocities at ^ that of the 
medium wave speed. The derivation details of the non-dimensional formulation are given 
in Appendix F. 

For experiments that follow, the reference wave speed is c 0 = 1, which allows the 
initial wave to propagate through the region of interest in a reasonable time for all exper¬ 
iments. Given the scale choices above yields a dispersion parameter f 2 of 0.1, however, 
to ensure that dispersion is felt, we choose f 2 = 0.5 corresponding to more than doubling 
of the angular velocity of the earth. Further, U and V are set to 0.0 or 0.1 under a two- 
dimensional Gaussian initial condition to test the formulation in a semi-infinite and infinite 
channel setting. 

The experiments begin with an analytic benchmark where a solution is synthesized 
to the KGE and compared to results obtained using the NRBCs. The chapter continues with 
more general experiments, which do not have analytic solutions. In channel experiments, 
two different initial conditions were tested. The first is a two-dimensional cosine pulse 
where b is the height of the channel (in all cases we used b = 4). This pulse is chosen with 
compact support such that the waves are generated only in a narrow region of the domain 
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and is zero elsewhere. Additionally, the initial condition (by design) satisfies the no-flux 
boundary conditions on the northern and southern boundaries. The initial condition is given 
by 

h(x, y, 0) = e~ 10x2 cos , h(x, y, 0) = 0. (VI.l) 

The other initial condition used in the channel and remaining experiments is a two-dimensional 
Gaussian centered at x = 0, y = 0, given by 

h(x, y, 0) = e - 10 (* 2 + 2 / 2 ) ) h(x, y, 0) = 0. (VI.2) 

In order to see the effect of the NRBC, we compare our solutions to the one com¬ 
puted on a larger domain (for all experiments except the analytic benchmark). We consider 
this reference solution “exact” when using equal order basis functions on a fine mesh, 
defined in greater detail in each of the subsequent experiments. Time integration was per¬ 
formed using various order Runge-Kutta schemes using a time-step chosen to ensure a 
Courant number of 0.25, where the Courant number is conservatively defined: 

Cn At 

Courant number = — ^ (VI.3) 

y (Ax) 2 + (Ay) 2 

Here, Ax and Ay are chosen as the minimum distance between any two points in the 
x or y directions respectively. This choice is made since the interpolation points are not 
uniformly distributed when using spectral elements. Consider a canonical element using 
8 th order basis functions with points chosen using the LGL interpolation points described 
in Chapter IV. As shown in Figure 12, these points are distributed in such a way to facilitate 
interpolation and integration and are not uniformly distributed. 


56 



+1 rr 


71 


_iU ;_: : : ttj 

-l +1 

Figure 12: Canonical element using 8 th order basis functions showing distribution of grid 
points. 


In order to quantify the errors observed between the reference and NRBC solutions, 
we use the normalized error defined as follows: 


error „ 

1 |ij n 



(VI.4) 


where N p is number of points in Q and and are the numerical and reference solutions 
at point k. What follows is a series of experiments that were designed to demonstrate the 
efficacy of the G-N NRBCs. 


A. ANALYTIC BENCHMARK SOLUTION OF SEMI-INFINITE HORIZON¬ 
TAL CHANNEL 


Recall the semi-infinite channel formulation discussed in Chapter IV that we now 
implement. Specifically, Y w is introduced at x\y = —2 and no-flux boundary conditions 
are specified on T^ and V s such that d y h = 0. Finally, the G-N NRBC B is introduced at 
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xe = 2 and advection in the x and y directions are set to zero. This specific situation is 
shown in Figure 13. 



Figure 13: The semi-infinite channel domain under consideration. Domain is truncated by 
an artificial boundary B at xe = 2. 


A priori, we synthesize a solution that satisfies the KGE with zero advection. Kucherov 
and Givoli [40] use a similar benchmark in a non-dispersive environment for a single wave 
mode. The solution used here is a linear combination of two waves and has the form 

2 

h bm (x, y,t) = cos cos (k m x - u m t + (j> m ) (VI.5) 

m =1 

where the parameters are the same as defined in (111.6). Given choices for k, b, c 0 , n, / and 
4>, we can determine c o to satisfy the KGE. Here, we choose /c = {|,|},6 = 4, Cq = 1, 
n = {2,4} , f 2 = 0.5 and 0 = {0,7r}. These choices ensure that the no-flux boundary 
conditions on T^ and Tg are upheld and the gradients on either side of the NRBC are 
matched at t — 0. For the NRBC solution, we specify the values for hr w based on the 
analytic solution. These parameter choices result in horizontal phase velocities Ck of 1.48 
and 4.22 for each of the waves. 

The NRBC parameters are selected simply as C 3 = c 0 , Vj as described in Chap¬ 
ter IV.D. Recall that this choice eliminates the second order time integration terms in the 
boundary formulation. In theory, if the Cj parameters were chosen to match the horizon- 
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tal phase velocities, J = 2 would be highest order NRBC required. In general problems, 
however, these phase velocities are not normally known a priori. We therefore keep our for¬ 
mulation general, relying on the fact that simply increasing the NRBC order is guaranteed 
to reduce the reflection caused by the boundary. 


1. Results 

In Figure 14, we show the reference solution on the top panel and the solution on 
the J = 4 NRBC truncated domain on the bottom panel at t = 3. The NRBC solution uses 
4 th order spectral elements on a 24 x 24 element grid, discretizing the domain into 9,409 
points. Qualitatively speaking, the results show very little reflection when compared to the 
synthesized solution. 

1.5 
1 

0.5 
0 

-0.5 
-1 
-1.5 


Synthesized Solution 



X 


NRBC Solution (J=4) 



X 


Figure 14: Semi-infinite channel comparing synthesized solution and the NRBC solution. 
4 th order spectral elements using NRBC order J = 4 with zero advection at t — 3 are 
shown. 


Quantitative results can be observed in Figure 15 showing the error on fl as a func¬ 
tion of spectral and NRBC order (J = 1,..., 10,15, 20). The number of elements is ad¬ 
justed for each spectral order to maintain an equal number of points (9,409) that the domain 
is discretized into. It is clear that increasing the NRBC order yields significant gains in ac¬ 
curacy, but by J = 5, the spatial discretization error dominates NRBC error in the low 
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order element (order 1 and 2) cases. However, by increasing the polynomial order of the 
elements, this spatial discretization error decreases and allows the true accuracy of the 
NRBC to be realized. 



Figure 15: Semi-infinite channel L'q error versus NRBC and spectral element order. Do¬ 
main is discretized into 9,409 points for all spectral element orders with zero advection at 

t = 3. 


B. SEMI-INFINITE HORIZONTAL CHANNEL 

For most problems, where initial data is generated from physical measurements, 
it would be impossible to generate an analytic solution with which to compare the NRBC 
solution. In this experiment, we use the initial data as described in (VI. 1) and (VI.2) to gen¬ 
erate the waves in the solution. To see the effect of the boundary condition, we compare 
our solution to one computed on a larger domain, i.e., — 2 < x < 10 where a homogeneous 
Dirichlet boundary condition h( 10, y, t) — 0 is prescribed for F E , replacing the NRBC. For 
this experiment, the discretization is chosen to maintain a mesh of 28,033 grid points for 
each spectral order. Time integration is performed with RK-8 to ensure the time discretiza¬ 
tion is not a limiting factor in computing the reference solution. We then solve the extended 
domain solution for t — 3, ensuring that the disturbance has time to propagate through the 
artificial boundary, yet has not had time to reach the eastern boundary. 
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Semi-Infinite Channel with Zero Advection Results 


In Figure 16(a), we plot the reference solution on the top panel and the solution of 
the truncated domain using the J = 4 G-N NRBC on the bottom panel. The solution on the 
truncated domain uses 4 th order spectral elements on a 24 x 24 element grid, discretizing 
the domain into 9,409 grid points. Qualitatively speaking, the results show very little re¬ 
flection using the combination of high-order elements and NRBC. Quantitative results can 


Expanded Domain Reference Solution 



NRBC Solution (J=4) 



(a) h(x,y, 3): U = 0, V = 0 


(b) £2 : u = o, V = 0 


Figure 16: Semi-infinite channel, 4 th order spectral elements (J = 4) using cosine pulse 
initial condition and zero advection. Left Plot: Contour plot showing h(x, y, 3) on ex¬ 
tended and truncated domains. Right Plot: Corresponding error versus NRBC and 
spectral element order. Domain is discretized into 9,409 points for all spectral element 
orders. 


be observed in Figure 16(b) showing the error on f) as a function of SE and NRBC order 
(J = 1,..., 10,15, 20). The number of elements is adjusted for each polynomial order to 
maintain an equal number of points (9,409) that the domain is discretized into. 

It is clear that increasing the NRBC order yields significant gains in accuracy, but 
by J = 5, the spatial discretization error dominates NRBC error in the low order element 
(order 1 and 2) cases. However, by increasing the spectral order of the elements, this spatial 
discretization error decreases and allows the true accuracy of the NRBC to be realized. 
Similar qualitative and quantitative results are found using the Gaussian initial data and are 
shown in Figure 17. 
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Expanded Domain Reference Solution 



(a) h(x, y, 3): U = 0, V = 0 (b) L 2 n : U = 0, V = 0 


Figure 17: Semi-infinite channel, 4 th order spectral elements (J = 4) using Gaussian initial 
condition with zero advection. Left Plot: Contour plot showing h(x,y,3) on extended 
and truncated domains. Right Plot: Corresponding Lq error versus NRBC and spectral 
element order. Domain is discretized into 9,409 points for all spectral element orders. 

2. Semi-Infinite Channel with Constant Advection Results 

In this section, we continue the analysis on the semi-infinite channel, this time 
adding constant advection in various directions. In this setting, we again must consider 
the selection of the Higdon parameters C r For the semi-infinite channel, we make our 
educated choice described in Chapters III.B and IV.D of augmenting the parameter value 
with the added advection. In the case of the eastern boundary NRBC, this yields a choice 
of Cj = c 0 + U which has the simplifying feature of making the boundary condition only 
first order in time. 

In Figure 18(a), we replicate the comparison plot between the reference solution 
and the truncated domain solution. Here, we use J = 4 G-N NRBC and 4 th order spectral 
elements on the same 24 x 24 element grid, this time with constant advection U = 0.1 
and V = 0. Qualitatively speaking, the results show very little reflection using the com¬ 
bination of high-order elements and NRBC. Further, if comparing this solution with the 
zero advection case, one notes the expansion of the wave in the direction of advection and 
compression where the wave is traveling against the direction of advection. These results 
match intuition and is the same behavior that van Joolen noted in [16, p. 108]. 


62 



















Expanded Domain Reference Solution 




(a) h(x,y, 3): U = 0.1, V = 0 


(b) L 2 n : U = 0.1, V = 0 


Figure 18: Semi-infinite channel, A th order spectral elements (J = 4) using cosine pulse 
initial condition with advection velocities U — 0.1 and V = 0. Left Plot: Contour plot 
showing h(x, y, 3) on extended and truncated domains. Right Plot: Corresponding L 2 Q 
error versus NRBC and spectral element order. Domain is discretized into 9,409 points for 
all spectral element orders. 


Quantitative results can be observed in Figure 18(b) showing the error on 0 as 
a function of SE and NRBC order (J = 1,..., 10,15, 20). The number of elements is 
again adjusted for each polynomial order to maintain an equal number of points (9,409) 
that the domain is discretized into. It is again clear that increasing the NRBC order yields 
significant gains in accuracy, but the spatial discretization error dominates NRBC error in 
the low order element (order 1 and 2) cases. 

Error is monotonically decreasing for each SE order, however, an interesting “dip” 
occurs for J = 3. Anomalies such as this can occur when choosing the C/s in a general 
manner as we have undertaken in this example. While the reflection coefficient guarantees 
that the reflection caused by the boundary will decrease as J increases, it says nothing about 
how much it will decrease. This depends heavily on the choice of the C/s; in this example, 
it is believed that our general C 3 choice for J = 3 happened to annihilate a significant 
wave mode or modes, resulting in better performance of the boundary condition. Similar 
qualitative and quantitative results are found using the Gaussian initial data and are shown 
in Figure 19. 
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Expanded Domain Reference Solution 



NRBC Solution (J=4) 




(a) h(x,y, 3): U = 0.1, V = 0 


(b) L 2 n : U = 0.1, V = 0 


Figure 19: Semi-infinite channel, 4 th order spectral elements (J = 4) using Gaussian initial 
condition with advection velocities U — 0.1 and V — 0. Left Plot: Contour plot showing 
h(x, y, 3) on extended and truncated domains. Right Plot: Corresponding L' A error versus 
NRBC and spectral element order. Domain is discretized into 9,409 points for all spectral 
element orders. 

As discussed in Chapter III, the Higdon NRBC implicitly assumes that any wave 
impinging on the NRBC is traveling primarily as a plane wave normal to the boundary. The 
previous example, where the advection velocity was in the same direction as the channel 
does not significantly challenge this assumption. In other words, to really test the value of 
the boundary condition, one must try advection velocities with some tangential component 
to the boundary. For these experiments, we consider only the Gaussian initial condition 8 . 

Figures 20(c) and 20(d) show the same L‘i A plots for advection velocities in other 
directions, one being the contrived case where the advection is perfectly tangential to the 
NRBC. These L' A plots correspond with the contour plots directly above that show the 
comparative solutions. Examining these results, it is clear that the boundary condition - 
even when put to a test with a wave pulse containing a significant tangential component 
to the boundary - performs well. It is noted that the order of the error suffers more under 
diagonal advection when compared to its individual axial counterparts. This may be due to 

8 Since the cosine pulse initial condition was constructed to ensure that the boundary condition oil I’ y and 
r s were automatically satisfied, any tangential advection would cause this condition to break down. For this 
reason, we choose to illustrate general results using only the Gaussian initial condition. 
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the additional terms activated in the interior and boundary formulations (IV. 9) and (IV. 10) 
when U and V are simultaneously nonzero. 

It should be further noted that as advection velocities are taken larger, the associated 
errors grow as well. In fact, experimentation shows that the formulation actually becomes 
unstable for high-order NRBCs as the velocities approach the reference wave speed c 0 . 
This is the same behavior noted by van Joolen in [31] when examined in a finite difference 
setting. This behavior does not seem to manifest itself, however, unless the advection 
velocities are taken much larger than would be expected in a geophysical setting. Further, 
this behavior seems to be exacerbated in all diagonal advection cases when using high- 
order (N = 16) basis functions. To counter this high-order instability, we implement high- 
frequency wave filtering as described by Giraldo et al. in [47, 50] when using N = 16 
basis functions. As in [47], we apply the filter every 10 time-steps at 20% strength. 

C. INFINITE HORIZONTAL CHANNEL 


For the next set of experiments, the domain is an infinite channel with the NRBCs 
located at x = ±2. The set-up is similar to that of the semi-infinite channel in that no¬ 
flux boundary conditions are specified on Tat and Y 5 . Further, T^ is exactly the same 
NRBC defined previously. This time, however, we prescribe another NRBC B on Yw- 
This specific situation is shown in Figure 21. 

The adjustments required for the western NRBC follow from the original derivation 
of the Higdon boundary condition as given in Chapter III.A. Now, instead of considering 
the right moving component of the wave approaching Ye, we consider the left moving 
component of the wave approaching Y w . To perfectly absorb the wave impinging on the 
boundary, we insist that F = 0 in (III. 10) such that the boundary satisfies the one-way 
advection equation h(x,t ) = G^x + (c 0 — U)t^j. The corresponding Higdon boundary 
condition is then given by: 


Hj : 




0 on r w ■ 


(VI. 6 ) 
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Expanded Domain Reference Solution Expanded Domain Reference Solution 



(a) h(x, y , 3): U = 0.0, V = 0.1 (b) h(x, y, 3): U = 0.1, V = 0.1 



(c) Lq: U = 0.0, V = 0.1 


(d) L 2 n \ U = 0.1, V = 0.1 


Figure 20: Semi-infinite channel, 4 th order spectral elements (J = 4) using Gaussian initial 
condition with advection velocities specified. Top Plots: Contour plots showing h(x, y, 3) 
on extended and truncated domains. Bottom Plots: Corresponding error versus NRBC 
and spectral element order. Domain is discretized into 9,409 points for all spectral element 
orders. 


In the experiments that follow, we use the initial data as described in (VI. 1) and 
(VI.2) to generate the waves in the solution. To see the effect of the boundary condi¬ 
tion, we compare our solution to one computed on a larger domain, i.e., — 6 < x < 6 
where homogeneous Dirichlet boundary conditions h(—6 , y,t) = 0 and h(6 , y,t) =0 are 
prescribed for Y w and T E , replacing the NRBCs. The discretization is again chosen to 
maintain a mesh of 28,033 grid points for each SE order. Time integration is performed 
with RK-8 to ensure the time discretization is not a limiting factor in computing the ref¬ 
erence solution. We then solve the extended domain solution for t — 3, ensuring that the 
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y N =2 


ys=-2 



Figure 21: The infinite channel domain under consideration. Domain is truncated by artifi¬ 
cial boundaries at xw = —2, and xe = 2. 


disturbance has time to propagate through the artificial boundaries, yet has not had time to 
reach the Dirichlet boundaries. 

1. Weak Form Adjustments 

Using a similar strategy for introducing a set of auxiliary variables for the western 
NRBC, applying them to the KGE equivalent, then converting any normal derivatives on 
the boundary to time and tangential boundaries, results in another, very similar formulation 
that is directly incorporated into the weak form (IV.3). The details of this formulation can 
be viewed in Appendix G. The selection of parameters Cj follows the “convenient” choice 
as developed for the eastern boundary, namely, to remove the second order in time auxiliary 
variable term by augmenting the reference wave speed with advection. This choice is 
Cj = Co — U for the western boundary. Similar experiments to those run in the semi¬ 
infinite channel were conducted using various advective velocities. 

2. Infinite Channel with Various Advection Velocities Results 

Qualitative and quantitative results using the cosine pulse initial condition are shown 
in Figure 22. In Figure 22(a), we plot the reference solution on the top panel and the solu¬ 
tion of the truncated domain using the J = 4 G-N NRBC on the bottom panel. Quantitative 
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results can be observed in Figure 22(b) showing the error on as a function of SE and 
NRBC order. In Figure 22(c), we replicate the comparison plot between the reference solu¬ 
tion and the truncated domain solution, this time with left to right advection. Quantitative 
results can be observed in Figure 18(b) showing the error on Q as a function of SE and 
NRBC order. In both cases, the number of elements is again adjusted for each polynomial 
order to maintain an equal number of points (9,409) that the domain is discretized into. 


Expanded Domain Reference Solution 



(a) h(x,y, 3): U = 0,E = 0 


(b) L 2 n : U = 0, V = 0 


Expanded Domain Reference Solution 





(c) h(x,y, 3): U = 0.1, V = 0 


(d) L 2 n : U = 0.1, V = 0 


Figure 22: Infinite channel, 4 th order spectral elements (J = 4) using cosine pulse initial 
condition with advection velocities specified. Left Plots: Contour plots showing h(x, y, 3) 
on extended and truncated domains. Right Plots: Corresponding L'i } error versus NRBC 
and spectral element order. Domain is discretized into 9,409 points for all spectral element 
orders. 
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Similar qualitative and quantitative results are found using the Gaussian initial data 
and are shown in Figure 23. It is noted that the order of the error again suffers more under 
diagonal advection when compared to its individual axial counterparts. 

D. OPEN DOMAIN CONSIDERATIONS 

In the construction of the G-N NRBCs presented thus far, we have assumed that all 
boundaries are aligned with the axial coordinates. Further, no two NRBCs have ever been 
placed adjacent to each other. In this section, we examine the consequences of placing 
NRBCs adjacent to each other in two configurations, namely the quarter plane and the 
open plane. This set-up is similar to the channel configurations described thus far in that 
the infinite domain is truncated via artificial boundaries B thus dividing the domain into a 
finite computational domain f) and a residual domain D. The only thing that changes is the 
configuration of the artificial boundaries. 

Specifically, the quarter plane is described as a domain that is bounded by physical 
boundaries T w and IV NRBCs are introduced at x = x E and y = y N . The physical 
boundaries are homogeneous Dirichlet conditions h = 0 on and T 5. This set-up is 
illustrated in Figure 24(a). The open plane is described as a domain that is unbounded on 
all sides. To compute a solution on such a domain, NRBCs are introduced at x = %, xe 
and y = ys, yN- This setup is illustrated in Figure 24(b). Artificial boundaries for F^ and 
T N are developed as outlined in Appendix G. 

1. Corner Compatibility Concerns 

To begin this discussion, consider the quarter plane. A source of concern is the 
method of handling the intersection point of T s and Tat. After all, the auxiliary variable 
form described in (III.20) and (G.9) are PDEs themselves and therefore require appropriate 
boundary conditions to be well-posed. In the channel, the no-flux conditions specified by 
the problem statement were applied to the auxiliary variables to make the problem well- 
posed. In the quarter plane, we have the homogeneous Dirichlet conditions for the western 
boundary of Tat and the south boundary of Te, but there are no such boundary conditions 
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Expanded Domain Reference Solution 


( 


□ 


NRBC Solution (J=4) 


(a) h(x,y, 3): U = 0,V = 0 


Expanded Domain Reference Solution 


LL 


NRBC Solution (1=4) 


(c) h(x, y, 3): U = 0.1, V = 0 


Expanded Domain Reference Solution 


tc—=■>! 


NRBC Solution (J=4) 


(e) h(x,y, 3): U = 0,V = 0.1 



(b) Li l :U = 0, V = 0 



(d) L z n : U = 0.1, V = 0 



(f) Lbu = 0,v = 0.1 



(g) h(x, y, 3): i7 = 0.1, V = 0.1 


(h) L z : U = 0.1, V = 0.1 


Figure 23: Infinite channel, 4 //l order spectral elements (J = 4) using Gaussian initial 
condition with advection velocities specified. Left Plots: Contour plots showing h(x, y, 3) 
on extended and truncated domains. Right Plots: Corresponding L q error versus NRBC 
and spectral element order. Domain is discretized into 9,409 points for all spectral element 
orders. „„ 
































































(a) Quarter Plane Domain 


(b) Open Plane Domain 


Figure 24: Left Plot: A semi-infinite domain Q truncated by artificial boundaries T E and 
Tat. Right Plot: An infinite domain Q truncated by artificial boundaries Ts, F^. T y and 

r w- 

that can be applied to the east of Tn and north of F/?- The question therefore arises, what 
are the appropriate boundary conditions for the auxiliary variables at these points? This 
problem is compounded in the open domain as there are no explicitly defined boundary 
conditions for any of the boundaries of the auxiliary variables. 

2. Use of Sommerfeld Radiation Boundary Conditions for Auxiliary Vari¬ 
able Boundary Conditions 

For this analysis, we suggest that the desired behavior of the boundary data on 
the auxiliary variables at these comers should minimize auxiliary variable reflection back 
into the computational (boundary) domain. In other words, the auxiliary variables should 
themselves be non-reflecting. Ultimately, we would like these boundary conditions to be 
easily implementable while still maintaining the true essence of the auxiliary variables. 
To implement this behavior , we consider a simple order J = 1 G-N NRBC (Sommerfeld 
condition) for the intersection points of two NRBCs, i.e., 

4>'Ax E ,y N ) = —— <j)j(x E ,yN) forT^ (VI.7) 

Co,</ 

1 . 

<t>'j{x E ,y N ) = - 4>j(x E ,y N ) for Tat (VL8) 

Co,a; 
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Here, c 0) y and c 0tX are the y and x components of the radial wave velocity and prime indi¬ 
cates the tangential derivative along the particular boundary, i.e., ^ along T E and 

(p' 0 = ^7 along Tat. We now recall the auxiliary form of the boundary condition for the G, 
terms onf E 


a j / Cifij-l dT E +Kj / Ciftj-l dT E — Ay / Cifij-l dFE + Pj / Ci4>j dT E 

Jr E dr E dv E dr E 

-7 [ Cifij dr E - f 2 f dr E = \ x f C 4 j+idr E . 

Jv E dr E dr E 

To ensure that the auxiliary form lies in H 1 (T E ), we integrate the second order in space 
term by parts to yield 


a,- 


dr E +Kj / dr E 7 x y / CiP'j-i dr E — \ y Q(f>'j 


Vn 


+ Pj / CiPj dri 


ys 


-7 / Ci0Wr B -/ 2 / c,o., 1 di' E ; a,. / GGj+i (ir, 


We now see that the boundary term contains (j)'- evaluated at the northern boundary. To 
implement the non-reflecting behavior, we make the substitution (VI.7) into the boundary 
term. The complete weak boundary form (for T E ) then takes the form 


T/ I QOj -1 c/T /?+kj / G0j—1 ^7'/■. 7 Ay / G0,—i ^F7 Ay GGj—1 7 0, / G Pj dT E 
d r E J r E J r E c o, y Jr E 

7 [ CiP'j dT E — f 2 [ Cipj-i dT E = X x [ (i(f>j + i dT E 


T E 


>Te 


T E 


for j = 1,_, J — 1; V G, 4>j € Vr B and <pj = 0 at ?/ s . A similar construct is readily 

computed for Tat. 

The only thing left to do is consider the problem of “double counting” the contri¬ 
bution at the comer. In other words, there are two values for the auxiliary variables at the 
comer; the one resulting from the evaluation of T E and the one from T N - Which aux¬ 
iliary variable contribution at the corner should be used, that of T E or that of T Y ? For 
this analysis, we adopt the “node-splitting” approach described by Pozrikidis [44, p. 215], 
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which amounts to averaging the corner auxiliary variable values. Of course, this formula¬ 
tion is inherently low order in the boundary treatment of the auxiliary variables. As such, 
we would not expect to see spectral convergence as shown in previous examples, however, 
there should be improvement over the J = 1 Sommerfeld condition. 

3. Corner and Open Domain with Zero Advection Results 

Figure 25 shows a series of contour plots showing how the initial disturbance prop¬ 
agates through the domain for 0 < t < 4.5 with zero advection. In this example, we run the 
simulation using 4 th order basis functions on a 24 x 24 - element grid (9,409 global points), 
using NRBC order J = 4. The simulation is run just long enough for the primary wave 
to exit the computational domain. Qualitatively speaking, the results appear to behave as 
desired - allowing the wave to propagate through the NRBC unimpeded. 

Time Evolution of the Quarter Plane Gaussian 


t=0.25 t=0.5 t=0.75 t=l t=l .25 t=1.5 



Figure 25: Time Evolution of quarter plane Gaussian (NRBC on T^ and T E ) using 4 Lh 
order spectral elements (J = 4) with zero advection. 
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Quantitatively, the results confirm that errors are not exponentially decaying as a 
function of NRBC order. As the order of the NRBC is increased, the crude approximation 
of the corner boundary condition on the auxiliary variables shows its weakness. In fact, 
experimentation shows that the boundary condition error quickly overtakes spatial and time 
discretization errors. Taking the Gaussian initial condition with tf = 3.0 for various basis 
function orders and boundary condition orders yields the L\ errors (using an extended 
domain solution as the reference) as shown in Table 2. 

Table 2: Lq Error as a function of NRBC Order for quarter plane using various spectral 
element orders. Gaussian initial condition and zero advection is used. 


NRBC Order 

Error 

Finear Elements 

L < } Error 
Order 2 Elements 

L'o Error 
Order 4 Elements 

J = 1 

0.11002 

0.10519 

0.10470 

J = 2 

0.07204 

0.04173 

0.05126 

J = 3 

0.04067 

0.03377 

0.03277 

J = 4 

0.03948 

0.02919 

0.03274 

J = 5 

0.03647 

0.02579 

0.02941 

J = 10 

0.03446 

0.02517 

0.02539 

J = 20 

0.03446 

0.02513 

0.02526 


Similar experiments were performed for the open plane. G-N boundary conditions 
are implemented along all four boundaries and the order 1 Sommerfeld boundary condition 
is applied to each auxiliary variable boundary as described in the previous section. Figure 
26 shows a series of contour plots showing how the initial disturbance propagates through 
the domain for 0 < t < 4.5 with zero advection. In this example, we run the simulation 
using 4 th order basis functions on a 24 x 24 - element grid (9,409 Global Points), using 
NRBC order J = 4. Again, qualitatively speaking, the results appear to behave as desired 
- allowing the wave to propagate through the NRBC unimpeded. 

Again, taking the Gaussian initial condition with tf — 3.0 for various basis function 
and boundary condition orders yields the L\ errors (using an extended domain solution as 
the reference) as shown in Table 3. Two main observations can be drawn from these results. 
First, the major source of error appears to be with the boundary treatment. While there are 
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Time Evolution of the Open Plane Gaussian 


t=0.25 t=0.5 t=0.75 t=l t= 1.25 t=1.5 



Figure 26: Time Evolution of open plane Gaussian (NRBC on all boundaries) using 4 th 
order spectral elements (J = 4) with zero advection. 

modest gains made by increasing the spectral element order, once J = 3, the errors for 
all spectral element orders are nearly the same. Second, even though the Liy results are 
less impressive than the channel experiments for high-order NRBCs, there is significant 
improvement from J = 1 (Sommerfeld condition) to higher J - even if the improvement 
is far from exponential. 

With this said, one may be concerned with the stability and long term behavior 
of this NRBC scheme employed in the quarter and open domain planes. Of course, as 
t —>■ oo, one would expect h 0. To gain a quantitative handle on this, consider the 
oo-norm defined as follows: 

IN loo = max \hi\ 
l<i< N p 

where N p is the number of points in Q. We choose this norm simply to get an estimate of 
how much of the initial disturbance is left in the computational domain after a substantial 
amount of time has passed. Using our now standard test case of 4 th order spectral elements 
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Table 3: Lf, Error as a function of NRBC Order for open plane using various spectral 
element orders. Gaussian initial condition and zero advection is used. 


NRBC Order 

IJf l Error 
Linear Elements 

Lq Error 
Order 2 Elements 

L’o Error 
Order 4 Elements 

J = 1 

0.23695 

0.21538 

0.21102 

J = 2 

0.09145 

0.08376 

0.04345 

J = 3 

0.06494 

0.04056 

0.04144 

J = 4 

0.06466 

0.03880 

0.03837 

J = 5 

0.06454 

0.03873 

0.03776 

J = 10 

0.06441 

0.03794 

0.03767 

J = 20 

0.06441 

0.03794 

0.03767 


on a 24 x 24 element grid with NRBC order J = 4, when computed for t = 1000 with 
J = 10, it was found that ||/z||oo = 1-03 x 10“ 17 forthe quarter plane and | |/i| |oo = 5.30 x 
HP 19 for the open domain; in both cases, essentially zero throughout the computational 
domain. While this is not a rigorous stability analysis, it does experimentally suggest a 
stable formulation. 

4. Corner and Open Plane Domain with Constant Advection Results 

It can be shown (see Appendix H for details) that when working in the open plane, to 
examine the behavior of the KGE under constant advective velocities in various directions, 
a simple change of coordinate system can recast the problem into a much simpler problem. 
The simplified problem is one where advection is in only one direction aligned with the 
new coordinate system. This implies that when examining the open plane under advection, 
it is sufficient to test only cases where advection is in the x or y direction. The benefits of 
this change of coordinate system include reducing the computational overhead, as well as 
minimizing various errors due to the more complex formulation if viewed in the original 
coordinate system. It should be noted, however, that the formulation discussed in VI.D.3 
still results in a stable formulation when diagonal advection is applied to the solution. To 
see this, examine Figure 27 where north-east advection is applied and the co, y and co, x terms 
have been adjusted by advection. 
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Time Evolution of the Quarter Plane Gaussian 

t=0.25 t=0.5 t=0.75 t=l t=l .25 t=1.5 



Figure 27: Time Evolution of open plane Gaussian (NRBC on all boundaries) using 4 ,/ ' 
order spectral elements (J = 4) with advection velocities U — 0.1 and V = 0.1. 


The corresponding L ^ errors are presented for various SE and NRBC orders in 
Table 4. 

E. EFFECTS OF HIGH-ORDER TIME INTEGRATION 

At the outset of this work, it was believed that at some point the improvements real¬ 
ized by increasing the spatial discretization and the order of the NRBC would eventually be 
limited by the time integration scheme [52]. To this end, the order of the time integration 
scheme was varied to examine the effects of time integration on accuracy of the solution. 
As has already been presented, gains made by increasing the order of the NRBC halt for 
lower order spectral elements after J = 5. Early experiments showed that even for high 
order (order 8 and 16) spectral elements, the gains made by increasing the order of the 
NRBC are limited at some point using classical RK-4 time integration. 
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Table 4: Lf, Error as a function of NRBC Order for open plane using various spectral 
element orders . Gaussian initial condition with advection velocities U — 0.1 and V = 0.1 
used. 


NRBC Order 

IJf l Error 
Linear Elements 

Lq Error 
Order 2 Elements 

L’o Error 
Order 4 Elements 

J = 1 

0.26914 

0.26529 

0.16998 

J = 2 

0.09571 

0.08642 

0.03810 

J = 3 

0.04391 

0.02711 

0.02631 

J = 4 

0.04061 

0.02514 

0.02467 

J = 5 

0.04028 

0.02413 

0.02460 

J = 10 

0.04021 

0.02410 

0.02460 

J = 20 

0.04018 

0.02409 

0.02460 


For this experiment, we consider the KGE on a semi-infinite channel with hp w = 0. 
To ensure that any boundary or time effects are not masked by the interior discretization, 
24 th order spectral elements are used on a fine mesh consisting of 4,753 global points. The 
Gaussian initial condition is used and is evaluated until t = 4. The reference solution in 
this case was computed as described previously, except this time using 24 th order spectral 
elements on a fine mesh consisting of 9,457 global points. Time integration was performed 
with a 10 th order Runge-Kutta scheme using a time-step chosen to ensure a Courant number 
of 0.1. 

As can be observed in Figure 28, gains made by improving the time integration 
matter only if combined with high-order treatment of the boundary. Conversely - gains 
using high-order treatment of the boundary can only be realized if there is a high-order 
treatment of the time integration. It should be noted that these results (error on the order of 
HP l0 ) cannot be observed unless high-order treatment of the interior also accompanies the 
high-order treatment of the boundary and time. Several experiments were conducted that 
varied the order of the interior, boundary and time integration [41, 53]. The clear result 
was that without high-order treatment of all components in concert, convergence to the 
reference solution is stalled. 
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Figure 28: Error in the SE Solution of KGE using NRBCs of various order as a function of 
time integration order. 

We believe as a practical matter that the order of all components of the numerical 
solution (spatial, boundary and time) should be chosen to ensure that the numerical method 
is the strongest (in accuracy) component of the model. If high-order treatment of any 
of the three components is missing, the high-order treatment of the other components is 
essentially wasted. For models where parameters and data have associated measurement 
and parameter errors, the numerical method should be chosen to maintain at least the same 
accuracy. 
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VII. TOWARDS ARBITRARY DOMAINS 


In the previous chapters, we examined the G-N boundary formulation using spectral 
elements on unstructured quadrilaterals. The physical boundaries, however were perfectly 
aligned with the coordinate system axes. This was convenient since it allowed the use of 
the G-N auxiliary formulation. The power of spectral elements lies not only in its ability to 
compute high-order accurate solutions, but also in its ability to handle complex geometries. 
While we demonstrated exponential error convergence in channel experiments when high- 
order treatment was applied to spatial, boundary and time components of the problem, this 
exponential convergence broke down when applied to boundary configurations where two 
NRBCs were adjacent to each other. In short, since there is a discontinuity in the normal at 
the intersection of adjacent NRBCs, the corner was the problem. 

This chapter considers what happens when we completely remove any corners. We 
first examine an arbitrary domain where the boundary can be of any shape. It will be 
shown that there are insurmountable complications that arise when using the G-N auxiliary 
formulation in this context. In this case, results are presented for various domains using a 
first order non-reflecting boundary condition and high-order G-N when certain simplifying 
assumptions are made. The chapter concludes by revisiting a boundary condition originally 
devised in 1998 for the wave equation by Hagstrom and Hariharan and modified in 2003 
by vanJoolen et al. to include dispersion. 

A. ARBITRARILY SHAPED BOUNDARIES 

Ideally, we would like to directly extend the work presented thus far to remove 
the problematic “corner” configuration and replace it with a continuous, smooth, closed 
boundary. If this were possible, then a single formulation (instead of four formulations 
combined in the open domain) for the boundary would result. The benefits to this type 
of boundary would be that the domain could be “fit” to the area of interest, reducing the 
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total number of grid points required. A general configuration that demonstrates this idea is 
shown in Figure 29. 



Figure 29: A general domain 0 truncated by artificial boundary F 


To begin, we again examine the Higdon boundary condition of order J as first pre¬ 
sented in Chapter III and the KGE as presented in Chapter II simplified by the assumption 
of zero advection. 


Hj: 




h 


h - CqV 2 /?. + f 2 h 


0 

0 


on r 


We note that the boundary condition and the KGE are described in two different coordinate 
systems - namely (n, r) and (x, y ) respectively where n and r are the normal and tangential 
directions on the boundary. If we consider an arbitrary part of the boundary (T) as shown 
in Figure 30, we can find a way to express the standard Cartesian derivatives in terms of 
normal and tangential derivatives. 

Of course, in the most general case, the normal and tangential vectors are depen¬ 
dent on the position on the boundary, i.e., n = ft(x,y ) and t = r(x,y). These normal 
components can be computed (see Appendix I) for a particular domain by considering a 
change of coordinates from (x, y) to (n, r) as defined by the linear transformation and its 


82 






n 



Figure 30: Components of normal and tangential derivatives 
associated differentiation operator: 
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Since the transformation is necessarily non-singular, this can also be written as: 
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(VII. 1) 


(VII.2) 


Now, if we expand the Higdon boundary condition, for J = 1, we get 


d n h + —h — 0 =>- n ■ V h = —— h. 

Li Li 


This is convenient since this boundary condition can be directly applied to the KGE zero 
advection weak integral form 


/ Vih dQ-c 2 0 / %fi ■ Vh dT + cl / W; ■ Vh dtt + f / %h dQ = 0 

In Jr Jn Jn 
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to yield the first order Sommerfeld formulation 


%h d<l - 

Li 


^ihdr + cl / Vtf* ■ V/i dfi +/ 2 / %hdn = 0. 


(VII. 3) 


/n 


1. Second Order (and Higher) Higdon Boundary Condition 

Recall that the Higdon boundary condition is very general. It can be applied to a 
variety of wave-type problems and reflection is guaranteed to decrease by simply increasing 
the order J. They suffer, however, from an implementation point of view since there are 
increasingly high-order spatial and temporal derivatives. Consider a second order Higdon 
boundary condition 

#2 : (d n + (d n h + -^-hj =0 on T. 

When expanded, this boundary condition takes the form 

/ 1 1 \ 1 .. 
d nn h + ( ~(j~ J (< (< 1 = 0- 


Continuing with the expansion to express the boundary condition in terms of the physical 
coordinate system, we find that d nn h is 


r\ r\ 

d nn h=—(d n h) = —(n-Vh) 
d ^ 

(fixOgh -|- rtydyh ) n ' V ( n x d x h riydyh ). 


(VII.4) 


The key point to take away from (VII.4) is that the components of the normal vectors are 
themselves functions of x and y. The product rule dictates that we must then compute 
the x and y derivatives of the normals in order to yield an “exact” representation of the 
higher order Higdon boundary condition. This expansion has two direct consequences that 
undermine efficient implementation of the formulation. 
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The first is something that has already been addressed, namely that the ability to 
implement the model is challenged - especially when dealing with large J. For each order 
that the Higdon formulation increases, high-order derivatives appear for h as well as the 
components of the normal vector. The second consequence is even more problematic in 
that the means to relate the boundary formulation back into the interior formulation has 
been lost. In the case of the Sommerfeld condition presented, the boundary condition 
and the boundary integral term (following integration by parts of the Laplacian operator) 
matched perfectly - thus allowing direct substitution of the boundary condition into the 
interior formulation. 

2. G-N on the Unstructured Boundary 

The G-N formulation was designed to remedy this problem of increasingly high- 
order derivatives by recasting them into a system of low order derivatives. If we try this 
with the unstructured boundary formulation, a similar auxiliary form to that presented in 
Chapter III.C) is obtained 

(d n + 0,-1 = 0,- j — (VII.5) 

where 0 O = h and 0j = 0. Again, the function h satisfies the KGE and 0i is obtained by 
applying the linear (constant coefficient) operator (d n + A-<9^ to h. 

Knowing that in the end, we would like to have an auxiliary variable formulation 
that contains only tangential and time derivatives (so that the boundary formulation can be 
discretized only on the boundary), we must consider the equation that 0, satisfies. It can be 
shown that when the KGE on the boundary is recast in terms of the normal and tangential 
coordinate system that it becomes a variable coefficient differential equation due to the 
presence of the normal components n x (x,y),n y (x,y). The result of this is that there is 
no general KGE-like equation that all 0, ’s will satisfy. In fact, every time we increase the 
order of the boundary condition, additional terms such as those encountered in (VII.4) are 
accumulated. In the end, a separate formulation that contains high-order derivatives will 
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have to be devised for each auxiliary variable introduced. In short, this implies that the 
G-N auxiliary variable formulation in its current form is incompatible with an unstructured 
boundary. 

To handle this incompatibility, we consider a simplifying assumption that the curva¬ 
ture on the boundary is small. This replicates the case where in a local region, the boundary 
“looks” like a straight line to the numerical solution. This assumption allows for all deriva¬ 
tives of normal components to be neglected, i.e., 


d_ 

dx 


(j^x) 



d . d . . 


The details of this boundary formulation and how it is integrated into the interior scheme 
are discussed in Appendix I. 


B. RESULTS FOR ADJUSTED G-N NRBCS ON ARBITRARY DOMAINS 


Some experiments were performed using the (what turned out to be) convergence 
limiting assumption of small curvature. While experiments showed stable behavior for the 
zero advection case, even over long term time integrations, the convergence was again, far 
from exponential in nature. For the first set of experiments, we consider how the formula¬ 
tion outlined in (VII.3) performs for various boundary shapes. Admittedly, this formulation 
is only a first order boundary condition, however as has already been shown, a first order 
condition is very easy to implement and has very modest computational overhead. The 
next set of experiments considers the adjusted G-N formulation. For these experiments, 
we consider rectangular, circular, and rounded rectangular boundaries where the Gaussian 
initial condition is used to generate the propagating waves. 

Figure 31 shows a series of contour plots for the zero advection case using 4 lh order 
spectral elements with J = 1 and J = 4. Model parameters are set to the standard test 
case where — l, f 2 — 0.5 and initial data as described in (VI.2) is used to generate the 
waves in the solution. To see the effect of the boundary condition, we compare our solution 
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to one computed on a larger domain, i.e., x,y e [—4,4] where a homogeneous Dirichlet 
boundary condition h — 0 is prescribed for I\ replacing the NRBC. For this experiment, 
the discretization is chosen to maintain a mesh of 28,033 grid points for each polynomial 
order. Time integration is performed with RK-8 to ensure the time discretization is not a 
limiting factor in computing the reference solution. We then solve the extended domain 
solution for t = 3, ensuring that the disturbance has time to propagate through the artificial 
boundary, yet has not had time to reach the boundary. The number of elements in each of 
the NRBC solutions is adjusted to ensure approximately 3, 000 global points were used. 

What is clear from these plots is that there are trade-offs between accurately rep¬ 
resenting the G-N NRBC and removing the problematic corners. Specifically, the square 
domain does not have to make an approximation for small curvature. In fact, with the 
exception of only 4 points (comers) in the global domain, the G-N NRBC is perfectly rep¬ 
resented by the arbitrary domain formulation. This is why there is significant improvement 
between the J = 1 and J = 4 cases. The rounded square and circular domains see less 
dramatic improvement as J is increased. While the problematic comers are removed in 
these cases, the small curvature assumption induces error in the G-N NRBC, thus imposing 
a convergence bound that cannot be overcome by simply increasing the order of the NRBC. 

A series of experiments was conducted to examine the normalized L'i A error as a 
function of spectral element and NRBC order for each of the NRBC boundary configu¬ 
rations. What is clear from the results shown in Tables 5-7 is that the errors are more a 
function of the NRBC than of the spectral element order. In short, there is almost no gain 
observed by using high-order spectral elements since much of the error is generated by the 
NRBC. When compared with the results shown in Chapter VI.D.3, which used the Som- 
merfeld approximation for the boundary condition of the auxiliary variables, the results are 
unimproved. 
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(a) Ref. Solution (b) Truncated Ref. Solution 



(c) NRBC J = 1 


(d) NRBC J = 1 


(e) NRBC J = 1 



(f) NRBC J = 4 


(g) NRBC J = 4 (h) NRBC J = 4 


Figure 31: Open Domain, A th order spectral elements using Gaussian initial condition with 
zero advection. Top Plots: Contour plots of reference solution solved on extended domain. 
Full and truncated domains shown for comparison. Center Plots: Contour plots of various 
NRBC boundary configurations using J = 1. Bottom Plots: Contour plots of various 
NRBC boundary configurations using J = 4 
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Table 5: Error as a function of NRBC Order for open plane arbitrary domain formu¬ 

lation using various spectral element orders on square NRBC domain. Gaussian initial 
condition and zero advection is used. 


NRBC Order 

Ij) Error 
Linear Elements 

Error 

Order 2 Elements 

L’o Error 
Order 4 Elements 

J = 1 

0.13632 

0.13475 

0.13561 

J = 2 

0.09030 

0.07383 

0.07164 

J = 3 

0.04483 

0.04403 

0.04293 

J = 4 

0.03996 

0.03993 

0.03841 

J = 5 

0.03964 

0.03984 

0.03841 

J = 10 

0.03948 

0.03957 

0.03841 

J = 20 

0.03948 

0.03952 

0.03840 


Table 6: L f } Error as a function of NRBC Order for open plane arbitrary domain formula¬ 
tion using various spectral element orders on rounded square NRBC domain. Gaussian 
initial condition and zero advection is used. 


NRBC Order 

L'q Error 
Linear Elements 

L\ Error 
Order 2 Elements 

L\ Error 
Order 4 Elements 

J = 1 

0.19805 

0.17845 

0.17331 

J = 2 

0.08245 

0.07027 

0.06646 

J = 3 

0.05528 

0.04839 

0.04458 

J = 4 

0.05246 

0.04599 

0.04358 

J = 5 

0.05194 

0.04547 

0.04234 

J = 10 

0.05187 

0.04540 

0.04203 

J = 20 

0.05186 

0.04540 

0.04201 


C. ALTERNATIVES 

Thus far, using the arbitrary boundary idea to remove the problematic comer points, 
we have not been able to improve results for the open domain problem. While the arbitrary 
boundary formulation does allow the user to choose a boundary domain of any shape (an 
advantage in certain contexts), the errors associated with this formulation were on par with 
alternatives presented in Chapter VI.D. Additionally, neither formulation led to exponential 
error convergence as the order of the NRBC was increased. With this in mind, we consider 
an alternative boundary formulation for a circular domain. Hagstrom and Hariharan [6] 
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Table 7: L\ } Error as a function of NRBC Order for open plane arbitrary domain formu¬ 
lation using various spectral element orders on circular NRBC domain. Gaussian initial 
condition and zero advection is used. 


NRBC Order 

L 2 ) Error 
Linear Elements 

I 2 j Error 
Order 2 Elements 

L’o Error 
Order 4 Elements 

J= 1 

0.26723 

0.26615 

0.25028 

J = 2 

0.10844 

0.10334 

0.09448 

J = 3 

0.07584 

0.07211 

0.06451 

J = 4 

0.07150 

0.06788 

0.06206 

J = 5 

0.07078 

0.06730 

0.06151 

J = 10 

0.07070 

0.06720 

0.06142 

J = 20 

0.07070 

0.06711 

0.06142 


devised high-order NRBCs for the standard time-dependent two-dimensional wave equa¬ 
tion in polar coordinates implemented in a finite difference setting. This NRBC follows 
the ideas pioneered by Bayliss and Turkel [1] with the exception that the NRBC condition 
does not involve any high-order derivatives after introducing auxiliary variables. Huan and 
Thompson implemented the same NRBC in a series of papers [54, 55, 56, 57, 58, 59] in 
a finite element setting. Here, we examine the effect of this work when examined with 
high-order spectral elements and time integration. 

1. Hagstrom Hariharan Polar Boundary Conditions 

The boundary condition devised by Hagstrom and Hariharan (hereafter referred as 
the HH formulation) provides a systematic approach for constructing boundary conditions 
for standard two-dimensional wave equation. The condition is based on the asymptotic 
series representation (which does not converge at any fixed radius) for an outgoing solution 
of the wave equation (in polar coordinates) 

1 d 2 h d 2 h 1 dh 1 d 2 h 

Cq dt 2 dr 2 r dr r 2 d6 2 ^ 

Since the boundary condition is asymptotic by nature, valid for large radial distances - this 
implies that larger radial distances should provide better NRBC convergence. Thompson 
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et al. make the observation that .. for practical problems, truncating the asymptotic ex¬ 
pansion after [ J] terms provides solutions with errors well below that of the discretization 
error” [59]. Here, we seek to significantly reduce the discretization error by employing 
spectral elements to find the true error convergence properties of the NRBC. In developing 
the boundary condition, Hagstrom and Hariharan construct a sequence of operators that 
approximately annihilate the residual of the preceding element in the sequence, viewed as 
a function on the artificial boundary. The sequence begins with a first-order Bayliss-Turkel 
operator discussed in [1], The boundary condition takes the form: 


dh 

dr 

4>j+i 


0i 


_1 _dh 

c 0 dt 


1 

2 r 


h, 


1 1 d^j-i 

c 0 dt r ' 3 4 r 2 J 1 4r 2 d6 2 




(VII.7) 
(VII. 8) 


where 


0o = 2 h and 0j = 0. 


At first glance, this boundary formulation suggests that we should develop a “new” 
spectral element formulation for the wave equation cast in polar coordinates. If we did this, 
however, we would then require a polar grid that would introduce additional complications 
such as the method of dealing with the degenerate quadrilaterals that inevitably occur at 
the center of the grid. Of course there are ways to overcome these obstacles, but it would 
be much more convenient to cast the problem in the same framework already developed. 
In other words, we seek to implement this boundary condition (presented in polar form) in 
our unstructured quadrilateral formulation of the wave equation (in Cartesian form). 

First, consider the two-dimensional wave equation (same formulation as presented 
in (11.29) with U = V = f = 0) 

S - cl^h = 0. (VII.9) 
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Multiplying by the test functions 'I', and integrating over the circular domain yields the 


weak integral form 


f) 2/ 7 


# ,;V 2 h (in = o 


Jn dt 2 Jn 

Transferring the second order spatial derivatives from h to the basis functions via integra¬ 
tion by parts and applying the divergence theorem to recast one surface integral term as a 
boundary integral gives us 


T JO 2 


Tkn • V h dn + Cn / VTq • V h dil = 0. 


(VII. 10) 


Of note now is that the boundary condition (VII.7) contains a radial derivatives of h that on 
the circle is precisely the normal derivative n ■ V/i. This allows direct implementation of 
the boundary condition into (VII. 10) as follows: 


r pdh r / 1 dh 1 \ f 

t'W dsl ~ c M + 4j™rVhdn = 0 . (vii. ii) 

Here, since on the boundary the radius is fixed, the ^ term may be treated as a constant. 

A similar weak form is constructed for the boundary formulation by multiplying 
(VII.8) by the test functions Q an d integrating over T yielding (after by integration by 
parts): 


I [ r d( p3 

co Jr ^ 9t 


dr + 3 - 

_i_^ d(f)j -1 end l r QQ 

4r 2 * 39 I start 4r 2 J r 86 d6 


dT 

dr 


Ci^j+i dr. 


(VII. 12) 


We now use the fact that the boundary is continuous and closed to surmise that the endpoint 
evaluation term vanishes. The formal problem statement is then: Find h G V and e Vr 
where j = 1,... J — 1, such that Equations (VII.11) and (VII.12) are satisfied V £ V 
and C; e W . 
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2 . 


Results for the HH Formulation 


A series of experiments was conducted to determine the effect of the HH boundary 
condition for various SE and NRBC orders. Since the formulation is designed for circular 
boundaries, we consider only circular boundaries with unstructured grids. In each case, we 
choose the number of elements to yield approximately 3, 000 global points. Since the Gaus¬ 
sian initial condition described in (VI.2) is perfectly symmetric with respect to the bound¬ 
ary, we introduce asymmetry by adjusting its shape to yield a smooth, two-dimensional 
“oval-shaped” initial condition with shape parameters a x = \,cr y = |, further rotated by 
an angle of 9 = |. The initial condition used here is: 

h(x, y, 0) = e-K+^W), h(x, y, 0) = 0. (VII. 13) 


Here, the parameters a, b, and c are defined as follows: 


cos 2 9 sin 2 6 

~ 2 a 2 + 2 a 2 

x y 


sin 29 sin 29 
4<r 2 4cr 2 

x y 


sin 2 9 cos 2 9 

2 a 2 + 2<r 2 

x y 


(VII. 14) 


Again, the solution is compared to one computed on a larger domain allowing the wave 
to propagate out of the NRBC domain but not yet impinge on the non-physical boundary 
used to compute the solution on the larger domain. Qualitative results are shown in Figure 
32 and quantitative L\ errors are shown for various NRBC orders for SE orders up to 6 in 
Table 8. No further improvement was observed for SE orders above order 6. 


3. Adjustments to HH to Include Mild Dispersion 

The unstructured grid representation of the HH formulation has been demonstrated 
to significantly reduce reflection caused by the boundary for the standard wave equation. 
The question now arises, can this formulation be extended to include dispersive effects such 
as Coriolis? In [60], van Joolen et al. presented a method to extend the HH formulation 
for the standard wave equation under mild dispersion. While this formulation was well 
grounded mathematically, as far as the author knows, it was never implemented. A brief 
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(a) Ref. Solution (t = 1) 


(b) Ref. Solution (t = 2) 


(c) Ref. Solution (t = 3) 



Figure 32: Open Domain, 4 th order spectral elements (J = 4) using oblique Gaussian 
initial condition shown for t = 1, 2, 3. Top Plots: Contour plots of reference solution 
solved on extended domain. Superimposed black circle indicates NRBC domain. Bottom 
Plots: Contour plots of various NRBC boundary configurations using J = 4. 

synopsis of their derivation follows with results presented for mild dispersion where f 2 = 

0 . 1 . 


We first consider the KGE without advection (in polar coordinates as in the HH 
derivation): 

1 d 2 h r) 2 h 1 dh 1 0 2 h f 2 

- J —h. (VII. 15) 


Cq dt 2 


dr 2 


1 dh 

r dr 


r 2 d9 2 Cn 


As has been previously discussed, in the geophysical context, the dispersion parameter is 
typically small. We assume here that 


c 0 K 


< 1 , 


(VII. 16) 
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Table 8: L 2 > Error as a function of NRBC Order for Hagstrom Hariharan NRBC formulation 
using various spectral element orders on the circular NRBC domain. Oblique Gaussian 
initial condition is used. 


NRBC Order 

L 2 ) Error 
Linear Elements 

L 2 ) Error 
Order 2 Elements 

L 2 , Error 
Order 4 Elements 

L 2 j Error 
Order 6 Elements 

J = 1 

0.09310 

0.04772 

0.04555 

0.04485 

J = 2 

0.03381 

0.00465 

0.00355 

0.00315 

J = 3 

0.02355 

0.00324 

0.00243 

0.00259 

J = 4 

0.02217 

0.00305 

0.00236 

0.00201 

J = 5 

0.02198 

0.00302 

0.00230 

0.00196 

J = 10 

0.02195 

0.00302 

0.00228 

0.00196 

J = 20 

0.02195 

0.00302 

0.00228 

0.00196 


where K is a typical wave number appearing in the solution. Now, apply the Fourier 
transform to (VII.6) and (VII. 15) in time to yield: 


cu 2 r d 2 h 


1 dh 
H— t; —f 


1 d 2 h 


= 0 


cu 2 f 2 \ ~ d 2 h 1 dh 1 d 2 h _ 

° dr 2 r dr + r 2 d6 2 


Wave 

Klein-Gordon 


where u is the frequency and h is the frequency domain representation of h. In both cases, 
we obtain the Helmholtz equation: 

k 2 h 

1 dr 2 ^ r dr r 2 dd 2 

In the non-dispersive case, K — — = K and K = \ K 2 — ^ in the dispersive case. In 

CO Y Cq 

order to facilitate the conversion back to the time domain, we now consider a Taylor series 
approximation to the square root term found in the dispersive case, i.e., 

X 

\Jl — x — 1 — -x + 0(x 2 ). 


95 




Provided that x is small, we can truncate the 0(x 2 ) terms. In our case, from (VII. 16) we 
can reasonably make this assumption yielding for the dispersive case: 


K I P~ f 2 - f 2 

_— , /1_ ■- _ ml - ■- _ K K _-_ 

K y CqA' 2 2c 2 IO 2 c 2 K' 

We now see that in the frequency domain, an equation valid in the non-dispersive case is 
valid in the dispersive case if we make the replacement: 


K ^\l K2 - f i~ K -2i- (VIU7) 

We now turn our attention to the boundary condition (VII.7) and (VII.8) that we Fourier 
transform in time to yield: 


—iK(j)j + -(f)j — 


-mil + 


dh 1 * 

d r 2 r 


{LlSl x _ 1 gjbzi 

4 r 2 f ^ 4 r 2 dO 2, 


01 

4 > j+ 1) 7 = Ij • • • 5 J 1 ■ 


(VII. 18) 
(VII. 19) 


Making the substitution (VII. 17), we obtain the dispersive version of the HH formulation 
in the frequency domain, i.e., 


- if' 2 j dh 1 - 
~ lKh+ 2M + fr + 2r h 


—iK(j)j + 


if 2 

2c 2 K 


V + “07 “ 


(./ )) : 

4r 2 


1 


1 <9 2 0j_i 
4r 2 <90 2 


01 

0j+i) 7 = !,•••,</ 1- 
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Transforming these equations back into the time domain results in the final HH 


boundary formulation for the KGE: 


i dh 

c 0 dt 


1 d(pj 

c 0 dt 


f 2 I' 1 i 

+ ^T / 0i(r)dr+-0 i 
zc o Jo r 

" -V-" 

n(t) 


f 


dh 


— / h(r)dT+ — + — h 


2 Ci 


o Jo 


U-i) 

4 r 2 


m(t ) 
2 


1 


<9r 2r 


1 

4r 2 


0i 

0j+l 


(VII. 20) 


(VII.21) 


where 


j = 1, • • •, J ~ 1, 0o = 2 h and 0j = 0. 


It should be noted that van Joolen et al. [60] show how m(t) and n(t) can be calculated 
in each time-step to keep the boundary condition local in time without having to store and 
operate on the history of the solution. For this analysis, a simple trapezoidal approximation 
was used to approximate the integral. 

The weak form of the formulation is now constructed. We consider the KGE in its 
general form: 

- CqV 2 /i + f-h = 0. 

Multiplying by the test functions 41,; and integrating over the circular domain yields the 
weak integral form 

/ (in - c 2 [ ^. t V 2 /i dtt + f 2 j V i hdn = 0 . 

Jn dt Jn Jq 


Transferring the second order spatial derivatives from h to the basis functions via integra¬ 
tion by parts and applying the divergence theorem to recast one surface integral term as a 
boundary integral gives us 


d 2 h 

* lW dn-cl 


4/jn • V/i dfl + Cq / V4q • V/i dCl + f 2 / dfl — 0. (VII.22) 
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Of note now is that the boundary condition (VII.20) contains a radial derivatives of h that 
on the circle is precisely the normal derivative n ■ V h. This allows direct implementation 
of the boundary condition into (VII.22) as follows: 




t Q2h 

V'W* 1 


Vi 


+ c o 


I f)h I f 2 \ 

*-«8i-* h -k* ) ) da 

Wi ■ Vh dQ + f I Vih dQ = 0. 

i Jn 


(VII.23) 


Here, since on the boundary the radius is fixed, the p- term may be treated as a constant. 

A similar weak form is constructed for the boundary formulation by multiplying 
(VII.21) by the test functions Q and integrating over T yielding (after by integration by 
parts): 


1 

c 0 


dp 


/ 


Ci-^dT + f- / Cin(t)dr+ J - / Ci^dT 


(J-IY 


dt 


2 c L ., 
4r 2 


<90j_i i end 
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r _ i 

1 

start 4 r 2 


^ 4 r 2 

<9Ci <90j-i 


<90 (90 


0T = 


Ci$j —i dr 
r 

Ci0j+i dr. 


(VII. 24) 


We now use the fact that the boundary is continuous and closed to surmise that the endpoint 
evaluation term vanishes. The formal problem statement is then: Find h G V and 0j e Vr 
where j = 1,... J — 1, such that Equations (VII.23) and (VII.24) are satisfied V Tq G V 
and 0 e V r . 


4. Results for HH with Dispersion 

A series of experiments was conducted to determine the effect of the HH boundary 
condition extended to include mild dispersion for various SE and NRBC orders. The set-up 
is identical to the experiments without dispersion, except the dispersion parameter is set to 
f 2 = 0.1. Qualitative results are shown in Figure 33 and quantitative Ip, errors are shown 
for various NRBC orders for SE orders up to 6 in Table 9. As in the non-dispersive case, 
no improvement was observed for SE orders above order 6. 
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(a) Ref. Solution (t = 1) 


(b) Ref. Solution (t = 2) 


(c) Ref. Solution (t = 3) 



Figure 33: Open Domain, 4 th order spectral elements (J = 4) using oblique Gaussian 
initial condition shown for t = 1, 2, 3 under dispersion f 2 = 0.1. Top Plots: Contour plots 
of reference solution solved on extended domain. Superimposed black circle indicates 
NRBC domain. Bottom Plots: Contour plots of various NRBC boundary configurations 
using J = 4. 
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Table 9: Error as a function of NRBC Order for Hagstrom Hariharan NRBC formulation 

using various spectral element orders on the circular NRBC domain. Oblique Gaussian 
initial condition is used with dispersion parameter set to / 2 = 0.1. 


NRBC Order 

L f, Error 
Linear Elements 

L't ) Error 
Order 2 Elements 

Lo Error 
Order 4 Elements 

L f } Error 
Order 6 Elements 

J = 1 

0.07290 

0.03555 

0.03369 

0.03293 

J = 2 

0.02684 

0.00371 

0.00283 

0.00248 

J = 3 

0.01869 

0.00258 

0.00192 

0.00204 

J = 4 

0.01759 

0.00243 

0.00186 

0.00157 

J = 5 

0.01744 

0.00240 

0.00181 

0.00154 

J = 10 

0.01742 

0.00240 

0.00180 

0.00153 

J = 20 

0.01742 

0.00240 

0.00180 

0.00153 
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VIII. CONCLUSIONS AND AREAS FOR FUTURE RESEARCH 


In this dissertation, we considered a reduced form of the shallow water equations in 
various (semi-infinite and infinite channels as well as open domain) configurations. Using 
the Givoli-Neta auxiliary variable formulation of the Higdon non-reflecting boundary con¬ 
ditions, we truncated the original infinite domain and developed the boundary conditions 
specific to the problem at hand. Using a high-order approach to the spatial discretization 
(spectral elements), time integration (high-order Runge-Kutta) in concert with high-order 
boundary treatment, we showed exponential convergence to the reference solution in chan¬ 
nel configurations. These results suggest a balanced approach to dealing with truncation 
errors - namely, to make improvements in all components of the problem to see improved 
accuracy. 

In open domain problems, we considered various ways to handle corner compatibil¬ 
ity concerns when using the Givoli-Neta auxiliary variable formulation. Using a physical 
argument that the auxiliary variables themselves should be non-reflecting at a boundary, 
we formulated a spectral element formulation that yielded stable solutions (even for long 
term time integrations) using first order NRBCs for the auxiliary variables. This formula¬ 
tion showed significant improvement from the first order (J — 1 Sommerfeld condition) 
to higher order J, although the improvement was far from exponential. Besides the low 
order method of handling the auxiliary variable boundaries, the “node-splitting” method of 
handling double-counting corner nodes turned out to be convergence limiting. 

Recognizing that any formulation that included corners would be problematic, we 
sought a boundary condition that would be valid for an arbitrarily shaped domain. The 
Givoli-Neta formulation was shown to have insurmountable implementation issues without 
the simplifying assumption of small curvature on the boundary. When the small curvature 
assumption was made, the formulation was shown to have stable, improved results from 
the first order (J = 1) Sommerfeld condition. It was clear, however, that there are trade¬ 
offs associated between accurately representing the Givoli-Neta formulation and removing 
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the problematic corners. Specifically, using a square domain (where the small curvature 
assumption affects only the four comer points) errors were shown to be improved over 
alternative domains (rounded square and circle). 

The final experiments conducted considered a boundary condition originally de¬ 
vised by Hagstrom and Hariharan and extended to the dispersive wave equation by van 
Joolen et al. This boundary condition is based on the asymptotic series representation 
of the wave equation in polar coordinates, valid for large radial distances. Results were 
improved over the alternative arbitrary domain formulations although valid only for large 
radial distances and restricted to circular boundary domains in this analysis. 

This research has demonstrated exponential convergence in channel experiments, 
only hypothesized in previous low-order settings. Additionally, it has developed several 
alternatives to handling open domains which improve performance to first-order NRBC 
schemes at a very moderate computational cost. What remains is to extend this high-order 
numerical formulation to more complex linear and non-linear systems of fluid motion such 
as the Euler equations. Additionally, better alternatives to dealing with corner compatibility 
concerns remains an open problem. 
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APPENDIX A. DEPTH INTEGRATING THE CONTINUITY 

EQUATION 


In order to arrive at the final form of the shallow water equations, we depth inte¬ 
grated our shallow water continuity equation. The details follow. Given 


0 = V ■ u 


we integrate in z 


0 


V • u dz 


' — hs 


f (t + t) 

J-h B \dx dy J 


dz + w\ z=h - w\ z= - hB 


(A.l) 


Since both h and h B depend on x and y (and t for h), we apply the Leibniz integral rule, 
which allows us to write: 


d_ 

dx 

d 


f‘Z=h 

rz=h 

du , 

dh 

/ udz = 


—— dz + u 

— - u 

/ z =—h B 

J z=—h B 

dx 

hdx 

nz=h 

rz=h 

dv 

dh 

/ vdz = 


-z—dz + v 

— - V 

J z=—h B 

J z=-h B 

dx 

hdx 


d{-h B ) 

h B dx 

d(-h B ) 

h B dx 


Substituting these into (A.l) we have 


d r h , 

0 = — / udz — 
dx J^h B 

d f h 

+ 7T / vdz - 

dy J- hB 


u 


dh 

z=h dx 

dh 

z=hdy 


— u 


d(-h B 
i=-h B dx 

d (— h B ) 

z=—h B dy 


+ w 


w 


z=h 


z=—hf 
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Now, using the appropriate boundary conditions 


w(x,y,-h B ) 
w(x, y, h, t) 


d h B d h B 

-U— - v- 


dx 


dy 


dh dh dh 

m +u di + l % 


we substitute to find 


d f h , d f h , 

0 = — / uaz + — / vdz 
9x J — h B dy J-h B 

dh 

+ di 


-u 


+u 


dh 

z=h dx 

dh 

z=h dx 


dh 
*=h dy 
dh 
t=hdy 


u 


u 


dh B 

z=-h B dx 

dh B 

z=-h B dx 


+ V 


dh B 

z=-h B dy 

dh B 
z=~h B dy 


Which simplifies to 


dh d f h d f h 

— + / udz + — / vdz. 

dt dx J_ hB dy J^ hB 


The construction of the shallow water model as shown in Figure 3 has H = h + h B where 
H is the depth of the fluid. Since a previous argument showed that u and v are independent 
of depth, we are left with our final, depth integrated continuity equation 


dh 

dt 


d_ 

dx 


(Hu) + 


d_ 

dy 


(Hv) = 0 
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APPENDIX B. LINEARIZING THE SHALLOW WATER 

EQUATIONS 


The non-linear version of the shallow water equations are: 

d t u + ud x u + vdyU — fv = —g d x h 
d t v + ud x v + vdyV + fu = —g d y h 
d t h + d x (Hu) + d y (. Hv) = 0. 

We wish to find a linear version of these equations. Suppose the bottom topography is flat 
such that h B is constant and u and v can be described by a constant mean term and a small 
0(5) deviation from that value, i.e., 

u = U + u* v — V + v* H — h B + h 

To be clear, U and V are the mean velocities and h B is the mean water depth. We now 
make these perturbation substitutions. 

dt(U + u*) + (U + u*)d x (u + u*) + (V + v*)dy(U + u*) - f(V + V*) 

d t (V + v*) + (U + u*)d x (V + v*) + (V + v*)dy(V + v*) + f(U + u*) 

d t h + d x (( h B + h) (U + u*)) + d y (( h B + h) (V + v*)) 

Now, recalling that U, V and h B are constants, we simplify to find: 

d t u* + (U + u*)d x u* + (V + v*)d y u* — f(V + v*) = —g d x h 

3 t v* + (U + u*)d x v* + (V + v*)d y v* + f(U + u*) = -g d y h 

d t h + h B (d x u* + d y v*) + d x h (U + u*) + d y h (V + v*) + h (d x u* + d y v*) = 0. 


= ~g d x h 

= -g dyh 
= 0 . 
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If we expand each term, the result is: 


d t u* + Ud x u* + u*d x u* + VdyU* + v*d y u* — f(V + v*) = —g d x h 
d t v* + Ud x v* + u*d x v* + Vd y v* + v*d y v* + f{U + u*) = -g d y h 
d t h + hsd x u* + hsdyV* + Ud x h + u*d x h + Vd y h + v*d y h + hd x u* + hd y v* = 0. 

We now neglect any terms of 0(d 2 ) to arrive at our final form of the linearized shallow 
water equations: 


d t u* + Ud x u* + VdyU* — f(V + v*) = —g d x h 
d t v* + Ud x v* + Vdyv* + f{U + u*) = -g d y h 
d t h + Ud x h + Vdyh + hs ( d x u* + d y v*) = 0. 
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APPENDIX C. ADJUSTING HIGDON’S CONDITION FOR 

ADVECTION 


Suppose that we have a wave that moves according to 

(d t + Ud x ) 2 h - c 2 0 dlh = 0. (C.l) 

This equation can be “factored” as follows: 

0 = (d t + (U — Co) dx'j (d t + (U + Co) d x ^jh 

Following a standard method of characteristics derivation, we define functions w and v as 

w(x, t ) = d t h + (U — c 0 )d x h 
v(x, t) = d t h + (U + c 0 )d x h. 

It can be verified that both 

d f w + (U — c 0 ) d x w = 0 
d t v + (U + c 0 ) d x v = 0 

satisfy (C.l) exactly. If we then solve w along its characteristic x = (U — c 0 ) t + x 0 and 

v along its characteristic x = (U + c 0 ) t + x 0 , with initial data w(x 0 , 0) = P(x 0 ) and 

v(x 0 , 0) = Q(x o), then the solutions for w and v are respectively 

w(x,t) = P^x — (U — co)tj = d t h + (U — c 0 )d x h (C.2) 

v(x, t) = q(x — {U + c 0 )t j = d t h + (U + c 0 )d x h (C.3) 
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We subtract (C.3) from (C.2) to find 

p(x — (U — c 0 )t^ = d t h P (U P c 0 )d x h 
- q(x- (U P c 0 )t ) = d t h P(U - c 0 )d x h 

P(x-{U - c 0 )tj -Q(x-(U + co)t) = 2c 0 d x h. (C.4) 

Further, if we combine (C.2) and (C.3) as 

(C - c 0 ) P(x -(U- c 0 )t) = (U- c 0 ) (dth P(UP c 0 )d x h ) 
— (U P Co) q(x — (U P co)f^ = (U P Co) (dth P ((/ — co)d x h^j 

(U - c 0 ) P (a - (U - c 0 )t) - {UP c 0 ) Q(x ~(U + c 0 )t) = -2 c 0 d t h. (C.5) 

This implies that the solution takes the form h(x, t) — F (x— ( U +c 0 )f j +G (x—(U—c 0 )t j. 
Here, F and G are arbitrary functions of the initial data. To see this, consider 

d t h = -{U P c 0 )F'(x -{UP c 0 )tj ~(U- c 0 )G'(x -{U- c 0 )tj (C.6) 

d x h = F' (x — {U P c 0 )f) + G' (x — (U - c 0 )t) . (C.7) 

Equating coefficients with (C.4) and (C.5) this yields a relation between the initial data and 
the functions F and G 

p(x — {U — c 0 )t) = 2c 0 G" (x-{U - c 0 )t) (C.8) 

q(x-{U P co)t) = -2 c 0 F' (x- (UP c 0 )t) . (C.9) 

This solution can be rewritten as 

h(x, t) = F (x - (c 0 PU)t) PG(xP (c 0 - U)tj (C. 10) 
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with the interpretation that the general solution is the sum of F yx — (co + U)tJ, a wave of 
fixed shape moving to the right with velocity (c 0 + U) and G (x + (c 0 — U)tj a wave of 
fixed shape moving to the left with velocity (c 0 — U). 
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APPENDIX D. NORMAL TO TANGENTIAL DERIVATIVE 
TRANSFORMATION FOR EASTERN NRBC 


The function h satisfies the dispersive, advective wave equation (11.29) in D. Since 
the function <pi is obtained by applying the linear (constant coefficient) operator [d x + ^rd t 
to h, it is can be shown that 0 \ should also satisfy the same equation in D 9 . Further, since 
(f)j is obtained by applying the same linear operator j — 1 times to (f> 1 , the functions <p :] 
should satisfy an equation like (11.29), namely, 

(dtt + ( U 2 — eg) d xx + (V" — eg) d yy + 

2 Ud xt + 2 Vdyt + 2 UVd xy + / 2 ) <j,j = 0 (D.l) 

Now, use the following identities: 

d xx 4>j = (d x — — dt j (d x + — dt j 4>j + -^2 

V L 'j +1 J V W+i / L y+i 

d x iOj Of ( 0 X Oj ) 

0 X y<Pj Oy ( d x (f)j ) 

and 

<9.x + J = -•*.> J - (D-5) 

9 Here we must use the assumption that Co and / are constants. By applying the differential operator to 
(11.29), computing each of the <f>j derivatives present in (III. 17) using the differential operator and simplifying, 
a simple induction argument shows that the (j>j’s must satisfy (III. 17) 


(D.2) 

(D.3) 

(D.4) 
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Now, if we substitute (D.2) - (D.4) into (D.l), and replace j with j — 1 everywhere 
yields, for j = 1,..., J 

<f>j -1 + (u 2 - Co) (ydx ~ -jj-dt'j $3 -1 + ^2 < A?- 1 + (^ 2 _ c o) dyy^j -1 + 

2 U dt (d x (f)j- 1) + 2V d y t4>j -1 + 2UVd y (d x <f)j- 1) + f 2 4>j ~i = 0 (D.6) 

From this and (D.5) one gets, for j = 1,..., J 

(f)j-i + ( U 2 - Co) ^ d x - 4>j + + (^ 2 ~ eg) dyy(t>j-i+ 

2Udt (d x (f)j- 1) + 2Vd y t(f)j-i + 2UVd y (d x (f)j- 1) + f 2 (j)j ~i = 0 (D.7) 

Now, we shift indices on (D.5) and multiply both sides by (f/ 2 — eg) as: 

(£/ 2 - eg) (d* + 4>i = (U 2 - eg) 0 J+1 j = 0,..., J - 1. (D.8) 

Subtract (D.7) from (D.8) 

1 + (^ 2 - c o) & ~ (^ 2 _ c o) ~ ( V2 ~ c o) dyvfa- 1 - 

2f/9 t (d*^) - 2Vd vt <!> j ^ - 2UVdy (c^-i) - /%_i = (^ 2 - eg) 0 i+1 (D.9) 

Now, consider (D.5). Expanding and solving for d x <j)j- 1 , we get: 

0*^-1 = (j)j - j — 1,..., J . (D.10) 
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Substitute (D.10) into (D.9) 


~ h -1 + (^ 2 “ c o ) ~ (^ 2 ~ c o ) 

(v 2 - Co) dyy(t>j -1 - 2Udt - 2Vd yt (t)j-i - (D. 11 ) 

2UVd y - f %.i = (f/ 2 - c 2 ) 0 J+1 

Simplify: 




U 2 — c o 


2C/1/ 


2K ^■_ 1 - (K 2 - c 2 ) 0"_ 1+ 


(V - Co) + gM - 2 uj i, - 21IV ^ - /Vj-, = (v 2 - cl) <f, j+1 (D.12) 




In (D.12) and elsewhere, a prime indicates differentiation with respect to y along Te, i.e. 
the tangential derivative along those boundaries. 


113 



THIS PAGE INTENTIONALLY LEFT BLANK 


114 



APPENDIX E. METRIC TERMS DERIVATION 


To facilitate interpolation and integration required for the SE method, we transform 
all terms of the weak integral form in physical space x = (x. y) 1 to a canonical space £ = 
(£, r]) T . This nonsingular mapping assumes x = x (£, rf) and conversely ^ = £ (x, y). 

A. DERIVATION OF METRIC TERMS 

Using the chain rule, we find 

, dx dx 
dx = — d£ + —dr/, 

which can be written in matrix form 


Here J e is the transformation Jacobian with associated determinant \J e \ defined as 


dx 

dx 

de, 

drj 

dy 

dy 

as. 

dr] 


dx dy dy dx 
<9£ drj d£ drj' 


Similarly, we can find the inverse transformation 


d£ = it dx + ^dy 
dx dy 


and write the derivatives of £(x, y) in matrix form as 
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where 


T e 


d£ d£ 

dx dy 
dy dr\ 
dx dy 


(E.4) 


is the Jacobian of the inverse transformation. Since the Jacobian described by (E.2) is 
non-singular, we can write the transformation described by (E.l) as 


dri 



Here ( J e ) 1 is the standard matrix inverse defined by 


(■jt 1 = 



dy _ ckc 

dy dy 

dy dx 


(E.5) 


(E.6) 


We now note that the formulations described by (E.3) and (E.5) are identical, therefore, 
equating coefficients from the two Jacobian terms (E.4) and (E.6) yields the metric terms 


<9£ 1 dy <9£ 1 dx dr] 1 dy dr\ 1 dx 

dx \J e \dr]' dy \J e \dy' dx \J e \ <9£’ dy \J e \ dC, 

This formulation is convenient since all metric terms are defined in terms of terms that are 
readily calculated via basis function expansions of cc(£, rf). 


B. CONSEQUENCES OF QUADRILATERAL GRID DEGENERATION 


Given the discussion of metric terms, we see that there is a common term in that has 
the potential to cause numerical instability - namely If the elemental Jacobian tends to 
zero, then all metric terms associated with that Jacobian will tend toward infinity. In fact, 
if the element Jacobian is small at certain points on the global domain compared with other 
locations on the global domain, experiments in this dissertation have shown that they tend 
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to corrupt the entire solution. The question then arises, what element geometries have the 
potential to cause numerical instabilities? 

Since the global problem is discretized into smaller elements, this question must be 
answered in the context of grid generation. There is general agreement in literature [45, 
61, 62, 63, 64, 65, 66] concerning quadrilateral grid generation that convex elements with 
maximum internal angles ~ 135° constitute a quality mesh. Li et al. in [67] describe pro¬ 
cedures to adjust quadrilateral basis functions to deal with elements that have large interior 
angles or, in fact completely degenerate into triangles, although with these adjustments, 
computational overhead is increased to deal with alternate canonical geometries. In this 
analysis, we seek a quality all quadrilateral mesh. 

To demonstrate how a poorly generated mesh can taint a solution, we consider an 
extreme example where a quadrilateral element is degenerated into a triangle as shown in 
Figure 34. In this case, one of the internal angles of the “ quadrilateral ” is 180°. 



Figure 34: Degenerate quadrilateral mapped to a canonical reference element. 
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Now, consider the linear basis function expansion (similar argument for higher or¬ 
der expansions) of the physical coordinates and their derivatives 

4 

Xn(£,v) = ^2^k(t, v) x k 

k= 1 

where the linear Lagrange basis functions are defined as 

V’i (V v) = I (1 - 0 (! - v), V >2 (V v) = | (1 + 0 (! - v) % 

^3 (V v) = \ (! - £) (! + v), A (V v) = \ (1 + 0 (! + v) ■ 

We consider the degenerate vertex located at (£, rj) = (1,1) and compute each term 
required for computation of the Jacobian. 

c)t 1 c)t 1 

^ (1, 1) = 2 (-*3 + * 4 ) , ^ (1, 1) = 2 (-*2 + * 4 ) , 

dv 1 dv 1 

^ (1, 1) = 2 (-1/3 + 1/4) , ^ (1, 1) = 2 (-V2 + y 4) . 

Using the fact that (x 2 ,2/2) , (^3,2/3) , (x 4 , y 4) are collinear, we put x 4 in terms of the 
other two points yielding 

*4= (y 4 -y3 + x 3 (^^-)) f^V (E.8) 

V \x 3 -x 2 JJ Vi/ 3 - 2 / 2 / 

Computing the determinant of the Jacobian (E.2) and simplifying using (E.8) yields 

\ J e \ = {- x 3 + X A ) (— 2/2 + 1/4) - (- 1/3 + 1/4) (~X 2 + X A ) ) = 0 (E.9) 

This degenerate quadrilateral will have infinite metric terms. Even if the quadrilateral el¬ 

ement does not completely degenerate into a triangle, but has a large angle - the metric 
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terms will be very large in comparison to other element metric terms and will have a desta¬ 
bilizing effect. Therefore, all meshes in this analysis were generated with internal angles 
less than 135° choosing to add additional elements (and degrees of freedom) to ensure that 
this happens. 
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APPENDIX F. NON-DIMENSIONALIZATION OF THE KGE 


Consider the KGE 


(d t + Ud x + Vdyf h - CgV 2 /r + fh = 0 


that we wish to study in a non-dimensional context. For simplicity in derivation, we assume 
that there is no advection (U — V — 0), yielding: 

S - + f-h = 0. 


Now, as outlined in [51], we examine typical scales of motion in the ocean so as to recast 
the problem in a dimensionless way (can substitute typical scales for atmosphere or other 
medium as well with the same process that follows). For this analysis, the length scales 
were chosen 0(100 km), vertical depth scales 0(100 m), scales for h 0(1 m) and the 
dispersion parameter / for Coriolis O (10“ 4 s _1 ). Given these choices, we know from the 
discussion in Chapter II that 

/ m \ m 2 

c n = 9 ^b = ( 10 — ) ( 100 m) = 1000 — 

V s 2 / s z 


that makes our specific problem: 


d 2 h 
dt 2 



V 2 h + 



0 . 


Now, we follow the details as outlined in [ 68 ] to scale out any dimensions. In particular, 
we define the following: 

x m 2 /m — h m 

x = — 7 — y = —7 — h =- 

10 5 m 10 5 m 1 m 
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where the typical length scale is 100 km (10 5 m). Now, we note via the chain rule that: 


dh 

dx 
d 2 h 
dx 2 


dh dx dh 1 
dx dx dx 10 5 m 
d_ (dh\ 1 d (dh 
dx V dx 


1 d (dh 


10 5 m<9o; \dx J 10 5 m9x \dx J 10 10 m 2 dx 2 


d 2 h 


similarly, 


d 2 h 


d 2 h 


dy 2 10 10 m 2 dy 2 


Using this information in (F) yields: 


d 2 h 1000 m 2 
dt 2 s 2 

d 2 h 


10 10 m 2 


(v 2 ft) 


10 ~ I „ 

H- ~^~h — 0 


1 v 2 h + ^h = o 


dt 2 10 7 s 2 

>10 7 s 2 ^ - V 2 h + 0.lh = 0 
dt 2 


Now, we remove the dimensions from our variable h to get 


10 7 s 2 ?-4 - V 2 h + 0.1h = 0. 
dt 2 


Letting t = |(| L ^ and noting via a similar argument as above that 


d 2 h _ 1 d 2 h 

dt 2 10 7 s 2 dt 2 


we arrive at our final, non-dimensional form of the Klein-Gordon Equation, 


d 2 h 

dt 2 


V 2 h + 0.1 h = 0 


where t = (10 3 -5 s) f, h = (1 m) h, x = (10° m) x and y = (10 5 m) y. 
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APPENDIX G. AUXILIARY VARIABLE LORMULATIONS LOR 
WESTERN, NORTHERN AND SOUTHERN BOUNDARIES 


For configurations studied in this dissertation, G-N auxiliary variable formulations 
are required for boundaries other than Te, explicitly derived in Chapter III. What follows 
here are the details of the formulation for each of the other boundaries. 

A. FORMULATION FOR THE WESTERN BOUNDARY 


We begin by stating the Higdon boundary condition for T w given by: 


Hj : 




0 on r iy . 


(G.l) 


When imposed as a boundary condition on Tw, we can recast this formulation in terms of 
auxiliary variables as outlined in Chapter III equivalent to the single boundary condition 
(G.l) as: 



where: © 0 = h 


j = (G.2) 


4>j = o 


This set of conditions involves only first-order derivatives. However, due to the appear¬ 
ance of the ^-derivative in (G.2), one cannot discretize the o 3 on the boundary Tw alone. 
Therefore, we shall manipulate (G.2) in order to get rid of the rc-derivative. 

The function h satisfies the dispersive, advective wave equation (11.29) in I). Since 
the function is obtained by applying the linear (constant coefficient) operator [d x — ^rd t 
to h, it is clear that 0 \ should also satisfy the same equation in I). Further, since 0 3 is ob¬ 
tained by applying the same linear operator j — 1 times to <p 1 , the functions o :/ should satisfy 
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an equation like (11.29), namely, 


(dtt + (U~ — Cq) d xx + {V~ — Cq) d yy + 

2 Ud xt + 2 Vdyt + 2 UVd xy + / 2 ) <j,j = 0 (G.3) 

Using the following identities: 

d X x4>j = (&x + 7/ dt \ (d x — — d t J (f)j + —^2 4>j 

V W+i / V W+i / L y+i 

dxt&j (c^.r 0j ) 

0/:y 07 0y ( 0r 0/ ) 


and combining with (G.2) allows us to write (G.3) as: 

(V - c 2 ) (4 + + 2(7) - 2UV<p'j - = (C/ 2 - c 2 ) * i+1 

for j = l,...,J-l (G.4) 

In (G.4) and elsewhere, a prime indicates differentiation with respect to y along T w , i.e., 
the tangential derivative along Tw As desired, the new boundary condition (G.4) does not 
involve ^-derivatives. In addition, there are no high-?/ or t derivatives beyond second order. 
The new formulation of the J th -order NRBC on T w can be summarized as follows: 

Poh + d x h = , 

ay0}-i + Kj$j_ i — + /3j0j — 70) — / 2 0j—i = 3>x0j+i (G.5) 

0o = h 0j = 0 
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where 


A) — 



a j 


2 U U 2 - eg 

q q 


Pi 



-2U, 


Kj = 


2UV 

~C~ 


— 2 V, 


7 = 2C/V, 


A, = K 2 - c 2 , 
\ X = U 2 - cl 


B. FORMULATION FOR THE NORTHERN BOUNDARY 


The Higdon boundary condition for T N is given by: 


Hj : 



(G.6) 


When imposed as a boundary condition on Tn, we can recast this formulation in terms of 
auxiliary variables as outlined in Chapter III equivalent to the single boundary condition 
(G.6) as: 




where: 


A/'-i 
</>o = h 


7 - (G.7) 


0J = 0 


This set of conditions involves only first-order derivatives. However, due to the appear¬ 
ance of the y-derivativc in (G.7), one cannot discretize the (p 3 on the boundary alone. 
Therefore we shall manipulate (G.7) in order to get rid of the y-dcrivativc. 

The function h satisfies the dispersive, advective wave equation (11.29) in D. Since 
the function (p 1 is obtained by applying the linear (constant coefficient) operator (d y + ^rd t 
to h, it is clear that <i>\ should also satisfy the same equation in D. Further, since <fij is ob¬ 
tained by applying the same linear operator j — 1 times to 7i. the functions o :i should satisfy 
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an equation like (G.3). Using the following identities: 


dyy(f)j ( Oy 


Ci 


3 +1 


9t I ( dy + 


C 


3 +1 


d t ) <j)j + 


3 ' 2 

0 i+i 


dyt (f)j {dy(f)j 


&xy*Pj 9 X (dy(pj ) 


and combining with (G.3) and (G.7) allows us to formulate the boundary as: 


(2V V 2 -c 2 0 ) 

{fa q J 


fa -1 + 


2 UV 

~fa 



- fa - (C / 2 - C 2 ) 4>U + 

- 2V ) 4 - 2f/V0'- - /Vj-i = 


(V 2 - c 2 ) 0 J +1 


for j = 1, ..., J - 1 (G.8) 


In (G.8) and elsewhere, a prime indicates differentiation with respect to x along T N , i.e., 
the tangential derivative along T \j. As desired, the new boundary condition (G.8) does not 
involve ^-derivatives. In addition, there are no high-rc or t derivatives beyond second order. 
The new formulation of the J th - order NRBC on T N can be summarized as follows: 


fah + d y h — <pi , 

fa—i fa— i T fa fa 70j / 4*j —i fa4 > j +1 (G.9) 

= h <j)j = 0 


where 


A) — 




G 2 


fa = (V 2 - c o) 



= 


2C/U 


-217, 


7 = 2 in/, 


A, = U 2 - c 2 , 

A X = U 2 - fa 
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c. 


FORMULATION FOR THE SOUTHERN BOUNDARY 


The Higdon boundary condition for Ts is given by: 






0 on r 5 . 


(G.10) 


When imposed as a boundary condition on T 5 , we can recast this formulation in terms of 
auxiliary variables as outlined in Chapter III equivalent to the single boundary condition 
(G. 10) as: 



C s ‘) 


where: 0 O = h 


j = !,..> ,-J (G.ll) 


<t>J = 0 


This set of conditions involves only first-order derivatives. However, due to the appear¬ 
ance of the ^-derivative in (G.ll), one cannot discretize the Oj on the boundary Ts alone. 
Therefore we shall manipulate (G.l 1) in order to get rid of the ^-derivative. 

The function h satisfies the dispersive, advective wave equation (11.29) in D. Since 
the function cf>i is obtained by applying the linear (constant coefficient) operator (0y ~ ±~d t 
to h, it is clear that 0 1 should also satisfy the same equation in D. Further, since <fij is ob¬ 
tained by applying the same linear operator j — 1 times to <t)\. the functions o :j should satisfy 
an equation like (G.3). Using the following identities: 

dyy(j)j 

dytrf’j 

^xy 0j 




Xj+\ 


dn — — -(9^ (j)j + 


c 


3 +1 


r 2 
u j +1 


dt (dy4>j) 

d x (dyCpj) 
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and combining with (G.3) and (G.ll) allows us to formulate the boundary as: 


( 2V V 2 -4 \ 
V Cj C) ) 



7-i - + 2 uj t'i-! ~ ( u 2 - eg) 4>U~ 

(0 + ch) + 2V ) X " 2UV4l i ~ /2 *"‘ = 7 2 " c ») ^+‘ 


for j = (G.12) 


In (G.12) and elsewhere, a prime indicates differentiation with respect to x along T s , i.e., 
the tangential derivative along IV As desired, the new boundary condition (G.12) does not 
involve ^-derivatives. In addition, there are no high-a; or t derivatives beyond second order. 
The new formulation of the J th -order NRBC onf s can be summarized as follows: 


l%h - d y h = , 

a j4>j -1 + K j$j -1 — A xfij-i + fij&j — 7 4"'j ~~ / 2( ^i-i = A y 0j+1 


3 

= h 


(G.13) 


0 j = 0 


where 


0 o — 

Pj = 


1 2V 

"Cx 


V 2 


C? 


= 


2UV 

~c~ 


2U, Xy = V 2 - cl 


( v 2 ~ (7 + 77) “ 2K 7 = 2UV ’ 


K = U — 4 
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APPENDIX H. OPEN PLANE DOMAIN ROTATION IN THE 
DIRECTION OF ADVECTION 


Here, we consider the effects of diagonal advection and how our formulation may 
be simplified. Suppose we have an open plane domain where NRBCs are specified on 
all four cardinal boundaries. The question arises - since we are dealing with an infinite 
domain, can the problem be simplified by a change in coordinates? If so, how would this 
adjust the problem, and would it make the problem any easier? 

To examine this question, recall our PDE in its standard x — y plane: 

(d t + Ud x + Vd y f h - CgV 2 /i + fh = 0 (H.l) 

Further, suppose that the advection velocities U and V are non-zero in both directions 
resulting in activation of all components of (H.l). To fix some ideas, let U >0 and V > 0. 
The goal is to convert (H.l) in its current coordinate system to one that is in the direction 
of the advection velocity. To see this, consider Figure 35. 


7/ y 



Figure 35: Generation of a new coordinate system in the direction of advection 
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The new coordinate system (£, rf) places the the £ axis in the direction of advection. 
This transformation is a simple rotation that can be described by: 

x = cos (0)£ — sm(9)rj 
y = sin(#)£ + cos(6)r] 


or: 


£ = cos (6)x + sin (0)y 
r] = — sin(6 ) )a; + cos {9)y 


Clearly U and V are also related by the geometry of the problem. 


sin(0) 

cos(0) 


V 

Vu 2 + v 2 
u 

Vu 2 + c 2 


(H.2) 


We note now that we can express h in terms of the new coordinate system as h(x, y, t) = 
h(£(x, y),r)(x, y),t). Since (H.l) is developed in the (x, y) system, we use the chain rule 
to expand (H.l) in terms of (£, rf). We adopt the shorthand convention 


_ dh _ 9,2 h 
''a "Tv 77T 

oa oaob 
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to yield the following expansions 


hx h v rix 

hy 7 Cy ~f 7 y 91/ 

hxii C^CrCy H - h^T] (JixTIy ~\~ ^y fjx) ~\~ hy 1 jT) x '9y 3“ ^l^Cxy "f" hyT) X y (H.3) 

hxx 7^4 — 9 O/s.r 7r H~ foyiiVx -f h^xx H - hyTjxx 

h%IV 7^4 A ^AyCydy -|- hyyTjy 4" ^l^yy A hyTjyy. 

Since the coordinate transformation is linear, the problem is simplified even more 
as the second order metric terms vanish. 

4 = cos 7 r r} x — — sin 0 £ y = sin0 r] y = cos(9) 

(H.4) 

Cra: Cyy Cry 9xx 9 yy 7/yy 0 

Now, looking term by term at (H.l), in light of (H.3), we consolidate terms to see 

the result of the transformation. 

7 + A h# + B7 W + Ch^ t + B/i^ t + E/i^ + f 2 h = 0 (H.5) 

where h = 7(£, 77 , t) and: 

A=(U 2 - cl) e x + ( v 2 - eg) C 2 + 2f/HC,4 
1 = (f / 2 - cl) f£ + (H 2 - eg) T/ 2 + 2UVr] x r]y 

C = 2 ([/£, + H4) (H. 6 ) 

B = 2 (Urjx + Vrjy) 

E = (f/ 2 - c 2 ) 2C x r/ x + (1 V 2 - eg) 24 ^ + 2(71/ (4% + £ v tj x ) 

We continue by using the information provided by the geometry of the transforma¬ 
tion in (H.2) and the metric terms found in (H.4). Specifically, if we define the adjusted 
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advection velocities U and V corresponding to the £ and // directions respectively as: 


U = U£ X + 

V = Urj x + Vriy 

then each of the terms in (H.6) simplify to: 

A = U 2 - C'l 
B = V 2 - C 2 0 

C = 2U (H.7) 

D = 2V 

E = 2UV - 2cq (2£ x r) x + i y r] y ) 

S "" ’v‘ ^ 

=0 


Of course, if one examines V we find: 


V = Urj x + Vrj y = —U sin(0) + V cos(6) = -U 


V 


+ ^- 


u 


Vu 2 + v 2 Vc 2 + v 2 


= 0 (H.8) 


This significantly simplifies the problem to: 


htt + {U~ — Cl ) h + 2 Uh^t + f 2 h — 0 


(H.9) 


This procedure shows that for the open plane domain, any constant advection ve¬ 
locity not in a cardinal direction in the standard x — y coordinate system can be converted to 
an equivalent problem where the advection velocity is in a cardinal direction in an alternate 
£ — rj coordinate system. The result is that when examining the open plane domain, one 
only needs to examine advections in one cardinal direction since a problem with diagonal 
advection could be recast into a cardinal direction advection problem in another coordinate 
system. 
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APPENDIX I. ARBITRARILY SHAPED BOUNDARY 
FORMULATION 


Recall the KGE 


(d t + Ud x + Vdyf h - c 2 0 V 2 h + f 2 h = 0 


(I-D 


for which we wish to formulate the G-N auxiliary variable formulation for an arbitrarily 
shaped boundary. The Higdon boundary condition of order J is given by 


Hj: 




h = 0 


on r 


( 1 . 2 ) 


Now, we introduce the auxiliary functions 0 1 ,..., 0j_ 1 > which are defined on T as well as 
in the exterior domain D (see Figure 29). Eventually, we shall use these functions only on 
T, but the derivation requires that they be defined in I) as well, or at least in a non-vanishing 
region adjacent to T. The functions 0. ; are defined via the relations 

i^dn + ^r h = 0i , (1.3) 

01 = 02 ) (1-4) 


(d n + dt ^ 0J_! = 0 . (1.5) 

By definition, these relations hold in D, and also on T. It is easy to see that (1.3 -1.5), when 
imposed as boundary conditions on T, are equivalent to the single boundary condition (1.2). 
If we also define 

0o = h 0j = 0 , (1.6) 
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then we can write (1.3 -1.5) concisely as 


(d n + 4>j- 1 = 4>j j = 1, • • •, J • ( L7 ) 

This set of conditions involves only first-order derivatives. However, due to the appearance 
of the normal derivative in (1.7), one cannot discretize the 4>j on the boundary T alone. 
Therefore, we shall manipulate (1.7) in order to get rid of the normal derivative d n = 
n x d x + riydy where n x = n x (x, y ) and n y = n y (x, y ). 

The function h satisfies the KGE in I). Since the function c>\ is obtained by applying 
the linear operator [d n + ^rdtj to h, we must consider what equation that 0 \ satisfies. We 
begin by applying this operator to the KGE to yield 

^<9 n C~ d ^j ^dxth + 2Vd yt h + 2UVd xy h + f 2 h^j = 0 

T X x d xxx h T X y d xyy h 2U d xx ^h 2Vd xy ^h -I- A x d xxy h H- f d x fij 

n y ipytth + A xdxxyh ^ydyyy h + 2 U dxyth + 2V dyyth + A x^xyy h + fdyh ) (1.8) 

7T ( dttt.h + A x d xxt h + A y d yyt h + 2Ud xtt h + 2Vd ytt h + A x d xyt h + f 2 d t h^j = 0 


Now, consider the first order derivatives for the <pi boundary condition 


<9*01 

<9y01 

<9 t 0i 


d_ 

dx 

d_ 

dy 

d_ 

Ft 


(n x ) d x u T fix&xxh + 

ijlx) (Xx'U + XlxfXxytl T 
( fi x ) d x u T n x d xt h T 


^ (^y) <9 yf-f fiyfX xy h ( t d x th 

d 1 

— yriy) dyU + Tlydyyh + —dyth 

d 1 

777 (ttj/) dyU -f 'YlyOyttl ()(f 11 

ot Cl 


Since the components of the normal vector are themselves functions of a: and y, we incur 
additional terms for the derivatives of those components. If the boundary is fixed in time, 
clearly the time derivatives of the normal components will vanish, however, the spatial 
derivatives will not. We now consider a simplifying assumption that the curvature on the 
boundary is negligible. In other words, the spatial rate of change of the components of the 
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normal vector are very small in comparison to the other terms. This allows us to neglect 
the following terms: 


dx M ~ dx^ y) ~ dy M ~ dy (ny) “ °’ 


Using this assumption, the second order derivatives of fa are 


& xx 4*i ^x^xxxh T H"yd X xyh T „ d xx fh 

t'i 

^xy 0 1 U/: dxxy ^ T Flyd X yyh 0 d X yt.h 

Or 

()yy0 1 U/: ()xyy h T~ ^y^yyy ^ T , , ^yyl ^ 

U i 

1 

(0t0 1 /r.y :Oxtth + flydytth 0 ^ ( dffifl 

Ul 

1 

9 x t.<P i r^ 7 ; 0 X xth 'f n y d xyt h T ., d x ^h 

v 1 

1 

dyt.tfil r/j; 0 X yi tl + Tlydyyfh T" 3'ytth 

Ui 


(1.9) 


We can now use (1.9) in (1.8) to find that C should also satisfy the KGE in I). Further, 
since (f)j is obtained by applying the same operator to 0j_i, the functions 4>j should satisfy 
a similar equation, namely, 


<9tt + (C 2 ~ Cq) + (V 2 — Co) d yy + 2 Ud x t + 2U d y t + 2UVd xy + /“ 


= 0 (1.10) 


We now note that the boundary condition and the PDE that 0j satisfies are described 
in two different coordinate systems - namely (n, r) and (x, y) respectively. Consider an 
arbitrary part of the boundary (T) shown in Figure 30. Of course, in the most general 
case, the normal and tangential vectors are dependent on the position on the boundary, but 
can be computed given a particular domain. Since these components then are “known,” 
we consider a change of coordinates from x and y to n and r as defined by the linear 
transformation and its associated differentiation operator (VII.2). Rewriting each operator 
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in (1.10) and again using the small curvature assumption yields: 


d tt = d tt 

&xx H'x^nn -tlj.H y() n r tlj f^ T r 
(hyy (lyd nrL “t“ — ^.rtty(hir 4“ ^x^tt 

9 x t Ylxdnt Tlyd T t 

dyt tty (h i t 4“ t t j . () T f 

thxy ti x tiy() nn 4" (tt x tty j d nT n x n y d TT . 

Substituting these terms back into (1.10) and organizing this simplifies to: 

[d u + Aa nn + Bd TT + Cdl + 2 m nt + 2 m Tt + f] y = o (i.ii) 

where: 

A = (U 2 - c 2 0 ) n 2 x + (V 2 - Cl) n 2 y + 2 UVn x n y 
B = (U 2 — cl) n 2 y + (V 2 - Cl) n 2 x - 2 UVn x n y 
C = 2 ( -U 2 + V 2 ) n x n y + 2UV (n 2 x - n 2 y ) 

D = Un x + Vn y 
E = — Un y + Vn x 

We now see that this formulation is expressed in terms of the normal and tangential coordi¬ 
nate system and further has the same form as (III. 17). We can proceed in the same manner 
as outlined in Chapter III to eliminate the normal derivatives to yield the new formulation 
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of the J th -order NRBC on F summarized as follows: 


M + §; 

ajh -1 + ij&j -1 ~ B 0J-1 + PAj - c 01 - f 2 4>j -1 




li - 7S 2 E > Pj - A 7T + 


3 

2D A 

“ ~ci ~ 1 “ 67?’ 


j 


C, 


cf a- 


-2D 


3 ^ 3 + 1 , 

00 = n 0J = 0. 


01 * 

A 0 J+ i 


( 1 . 12 ) 

(1.13) 


(1.14) 


Here, prime indicates tangential differentiation along T. As desired, the new boundary 
condition does not involve any normal derivatives and there are no high tangential or time 
derivatives beyond second order. 

There are a few remaining concerns. First, we must be able to compute the normal 
and tangential components of the vectors mandated by the mapping (VII. 2). Addition¬ 
ally, we have several terms which require the integration of tangential derivatives along a 
general boundary. The question arises, how do we evaluate these terms along a particular 
boundary? Finally, we must still relate this boundary formulation back into the interior 
formulation and consider the appropriate values for the C 3 terms. 

A. COMPUTING THE NORMAL AND TANGENTIAL VECTOR COMPONENTS 

Consider (VII.2), which gives us a way to write the normal and tangential deriva¬ 
tives in terms of the standard x — y coordinates. We can extend this via the chain rule to 
map each of these components in terms of our canonical £ — r/ coordinates. Before we do 
any of this, however, we must first be able to compute the normal vectors at each point 
on the boundary. To see this, consider the normal and tangential vectors of a canonical 
element as shown in Figure 36. 
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Figure 36: Normal and tangential vectors along a canonical element boundary. 

After introducing the nonsingular mapping £ = £(x, y) and rj = rj(x,y), these 
normal vectors [47] are given by: 


n = 77 


V /7 


|Vjj| 


r)=± 1 


n = £ 


V£ 


|V£| 


?=±i 


Each of these normal vectors can then be used to compute the tangential vectors by 
taking the cross product of the normal vectors with the unit vector (0,0, — 1) T . In the case 
of Figure 36, the normalized tangential vectors are: 




When put in terms of (IV. 8 ), the unit tangential vectors are: 


Ti = 



t 3 = 





(1.15) 
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where: 



We now have a way to compute the normal and tangential components required for this 
formulation. 

B. INTEGRATION OF TANGENTIAL DERIVATIVES 


In (1.13), we have several terms that will require the integration of first and second 
order tangential derivatives along a general boundary (after weak integral formulation). 
The question arises, how do we evaluate these terms along a particular boundary? 


1. Integration of First Order Tangential Derivatives 

To examine this arbitrary domain formulation in greater detail, consider a single 
element where we evaluate an integral that contains a first order tangential derivative along 
a single canonical side. 



d T = f- V 


Q = number of quadrature points on a side 


Here, J*| is the Jacobian of the transformation along the side (side length) and w q is the 
quadrature weight for Lobatto integration. If we use collocated quadrature, the cardinality 
property of the basis functions ensures that only the i th basis function is nonzero (in fact 
unity) at quadrature point i. This rids us of the summation over all quadrature points for 
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basis function i yielding: 


, d(p{x,y) | TS , ( i d(p i , idft 

iT ' = TO,|Ji 1 v* &T + T » aj" 


Now, expand our variable 0 in terms of our canonical coordinates (£, //) and perform basis 
function expansions using our 2-D basis functions (Mjv = (N + l) 2 ). 


Mjv 


Mjv 






j=i 

_ ^ <90* <9£* + dtp] dr) 


. 7=1 ^ 


'x / , 
J =1 


<9£ <9a: dq dx 


Mjv 


+ t jE 

i=i 


dtp] dp 1 dtp] drf 


dp dy dr) dy 


Of course our boundary integral implies that we are integrating strictly on a side. For 
definiteness in the example, consider side 2 where p = +l,i) £ [—1, +1]. Now, we expand 
using (IV.8) and (1.15): 



Simplifying and using the definition of the element Jacobian (IV.7), our task now requires 
evaluating: 


=Wi 



/ dx l dy 1 dy 1 dx 1 \ yX, dtp] \ 

\ dr) dp drj dp ) dr) 3 J 


Mjv 


= Wj 


Q r , 

.7 = 1 


Qipp 

Now, consider the term -0, which means “the partial derivative with respect to // 
of the j th basis function evaluated at the i th quadrature point”. We note that the 2-D basis 
functions tpj are formed using tensor products of 1-D basis functions ly-ip) and vi(rj) such 
that j = 1,..., M n where the mapping from 1-D to 2-D is j = {(k — l)(iV + 1) + l : 
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k, l — 1,... N + 1}. Therefore: 




Using this and the cardinality property of the basis functions, we can evaluate the integral 


of the first order tangential derivative as: 


f , d<!>{x, y ) ^ du ) 


A similar argument holds for other canonical sides being evaluated on the boundaries. 

2. Integration of Second Order Tangential Derivatives 

Recall, the second order term from the auxiliary formulation of the form: 


/ JT f 9 f d( t ) \ JT f jt 

dr - =X, §? J dr * - X, * 

dd end r d ty. d(f) 

f* dr start y rs <9r <9r d s ' 

=0 if closed boundary 


Now, suppose that we expand each term, as in the previous section. 


/• ch/y deb f ( dipt dx/ji\ ( <90 d<p\ 

L dr * = X. U & + T ” % J Xs + T *Ty) dT ‘ 


Q 

q= 1 




•s| r q i 

9 V s dx y dy 


d(j) q d(j) q 
Tx lkc +Ty l)y 
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Consider ® evaluated on side 2. 


r T y 


dx 


dy 


^( d^jdj 9 dxjil drf \ q / d£ q dh/{ drf \ 

x \ <9£ 9a; dy dx J Ty \ <9£ dy dr) dy ) 

1 dx f d'df ( 1 dy \ d'dl f 1 <9y \ \ 

\Jq\dr) \ <9£ VI® I drj) + dy V \ J q\ d UJ 
1 dy f dt{/- / 1 dx\ d'df / 1 <9a;\\ 

Jq\dy \ \je\dr,) + dy Vl J g e l^// 

1 / dxdy dy dx \ dh/- 

\Jq\\ Jq I V dy <9£ + dy d£) dy 
1 

I Jq I dy 


Now, combining this result with the expansion of (2) as found in Appendix I.B, we find: 



As in Appendix I.B, similar results can be obtained for each of the other canonical sides. 


C. RELATING THE BOUNDARY AND INTERIOR FORMULATIONS 


The final four terms of (IV.3) are all boundary integrals where the NRBC must be 
applied. To simplify our discussion, we note that these integrals only apply to h on T, thus 
we will denote these terms as h. Much like the auxiliary variable formulation, we note that 
the NRBC is defined in terms of normal derivatives (d n ) while the boundary terms of (IV.3) 
are defined in terms of standard Cartesian derivatives (d x and d y ). Again, we consider a 
linear transformation of the final four terms as defined in (VII.2). This yields, for the first 
boundary term: 

n x dF = \ x J T* 


dh f dh 

^ixr~n x dT = \ x / r,— 
OX Jr OX 


dh dh 

n x — - n.„— I n x cii (1.16) 


' dn 


'dr 
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Now, we can directly use (1.3) to join the two formulations: 


A 


X 


% 


dh df) \ 

x dn Uy dr J 


n x dT 





n x dT (1.17) 


Now, we note that the boundary integral can be discretized on the boundary alone - in 
terms of tangential derivatives, time derivatives, and the auxiliary variable di. Performing 
these substitutions in each boundary term of (IV.3), we get the revised weak form of the 
problem: 


r ji dQ~ X x I 

Jo OX UX 


d% dh da - a. 


d^i dh 

d\L 


dy dy 


-UV 


O'!/ dh 

1 drt-uv 


dy dx 


d^i dh 

d\L 


dx dy 


dh 


dh 


+2U / r t —dn + 2v / %—dn + r / %hdn 


' dx 


dy 


dh 

4 -A^ / r.idin 2 x dr - X x j f3 0 n 2 x dr - X x j 'ff i —n y n x dr 


(1.18) 


+A y j ^ i4>in 2 y dr - Xy J rji (3 0 n 2 dr + Xy J r~n x n y dr 
+UV Jr idin x n y dr — UV J rih/3 0 n x n y dr — UV J ^iJJti 2 dV 

+UV I r i<j)in x n y dr — UV I rih, /3 0 n x n y dr + uv f dr = 0, 

Jr Jr Jr dr 


which, as desired, evaluates boundary integrals with data derived from only the boundary. 


D. SELECTION OF Cj TERMS 


If we examine the auxiliary variable formulation (1.13), we see that the selection of 
appropriate Cj values has yet to be addressed. As has been previously shown, any choice 
of Cj is guaranteed to reduce spurious reflection as the order of the NRBC ( J) increases. 
Armed with this, we choose convenient values for our Cj ’s that cause the second order in 
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time (ay) terms to vanish for all j. In this case: 


„ 2D 1 A 

aj ~ 0 “ a ~ 1 “ c? 


=> C, = D ± vD - A 


= Un x + Vn y ± ^eg (n2 + n*) 
= Un x + ± Co 


This implies a special value of Cj for each boundary point, independent of NRBC 
order, i.e., C 3 = C(x. y) for all j. Further, this implies that the other coefficients in the 
auxiliary formulation are independent of NRBC order: 


C 

C(x,y) 


-2E, 


Pj = P 


2A 

C(x,y) 


-2D 


The Higdon boundary condition, while general in nature, implicitly assumes that by the 
time a wave pulse gets to the artificial boundary, it is traveling primarily as a plane wave 
normal to the boundary. This choice for C 3 can be thought of as a choice that accounts for 
any advection and “corrects” for the geometry of the problem - i.e., non-normal impinge¬ 
ment on the boundary. 
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